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SUMMARY  ..  _ _ 

During  the  year  of  1986-1987 ,  we  have  developed  the  vortex  element 
method  and  the  transport  element  method,  for  the  numerical  simulation  of  the 
NavierAStokes  equations  and  the  energy  and  species  conservation  equations, 
respectively.  These  methods  are  based  on  a  Lagrangian,  grid-free,  time- 
accurate  simulation  of  the  governing  equations  at  high  Reynolds  and  Peclet 
numbers,  without  resorting  to  turbulence  modelling.  Finite^tate  chemical 
reactions,  finite  compressibility  and  finite  heat  release  rates  are  also 
considered  in  the  formulations  of  the  numerical  schemes.  To  validate  these 
methods,  we  are  obtaining  solutions  for  reacting  shear  layers,  both 
homogeneous  and  heterogeneous ,  under  various  idealizations,  and  comparing 
the  numerical  results  with  experimental  data.  The  solutions  are  also 
analyzed  to  investigate  the  mechanisms  of  turbulence^combustion 
interactions.  ,  — 
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OBJECTIVES 


The  goal  of  this  research  program  is  to  develop,  test  and  apply  methods 
of  numerical  simulation,  based  on  vortex  dynamics,  to  reacting  flows  with 
finite  chemical  reaction  and  heat  release  rates.  In  particular: 

I.  The  development  of  accurate  and  efficient  numerical  methods  which 
can  be  utilized  in  the  simulation  of  the  time-dependent,  multi¬ 
dimensional  Navier-Stokes  equations  at  high  Reynolds  number,  and 
can  be  extended  to  solve  the  energy  and  species  conservation 
equations  in  cases  where  the  chemical  reaction  rates  are  finite  and 
fast  and  when  the  associated  heat  release  is  large  and  hence 
dynamically  important. 

II.  The  application  of  these  numerical  algorithms  to  turbulent  reacting 
shear  layers,  for  both  homogeneous  and  heterogeneous  reactants,  to 
validate  the  numerical  methods  against  experimental  results,  and  to 
study  the  underlying  mechanisms  of  entrainment  and  mixing  and  how 
they  affect  the  rates  of  product  formation. 

III.  The  investigation  of  the  mechanisms  of  turbulence-combustion 
interactions  based  on  rigorous  fundamental  models,  and  how  these 
interactions  can  be  manipulated  to  provide  more  efficient  burning. 
In  particular,  the  effect  of  turbulent  fluctuations  and  flow 
stretch  on  the  rate  of  chemical  reaction,  flame  stability  and 
extinction. 
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Three  graduate  students  have  completed  their  master's  theses  under  the 
sponsorship  of  this  program.  Their  names  and  thesis  titles  are  listed  at 
the  end  of  the  report.  Currently,  five  graduate  students  are  continuing 
their  doctorate  work  under  the  partial  of  full  support  of  this  project. 
Their  names,  listed  according  to  seniority,  are 

(1)  Ghassem  Heidarinejad 

(2)  Omar  Knio 

( 3 )  Habib  Na jm 

(4)  Anantha  Krishnan 

(5)  Luis-Fillipe  Martins 

EQUIPMENT 

To  meet  the  computational  needs  of  this  work,  we  have  built  the 
following  system  around  a  VAX  11/750: 

(1)  An  array  processor  MAP-6420. 

(2)  Two  MicroVAX  II  workstations,  and  a  MicroVAX  cpu. 

(3)  A  local  area  network  with  communication  interface  to  a  supercomputer. 

(4)  Graphics  terminals,  a  laser  printer,  and  a  color  film  recorder. 

The  system  has  been  hard-wired  into  the  Campus-wide  network  to  allow 
easier  access  to  other  computational  facilities  available  within  and  outside 
M.I.T.,  especially  the  John  von  Neumann  supercomputer  center  at  Princeton. 


In  the  following,  the  vortex  and  transport  element  methods  are  briefly 


described  and  the  results  of  their  preliminary  applications  are  discussed. 
THE  VORTEX  ELEMENT  METHOD 

In  this  method,  Lagrangian  particles  are  treated  as  finite  vortex 
elements  that  accurately  discretize  the  vorticity  field.  We  have  shown  that 
the  accuracy  of  the  method  is  governed  by:  the  shape  function  of  individual 
elements,  the  core  radius  and  the  distance  between  neighboring  elements.  To 
preserve  the  accuracy  as  the  flow  develops  strong  strain  fields,  particle 
distribution  must  change  to  accommodate  distortions  of  the  vorticity  by  the 
strain  field,  leading  to  a  natural  growth  in  the  number  of  vortex  elements 
with  time  as  the  flow  develops  stronger  gradients,  or  fine  scales,  via 
stretch.  We  are  implementing  a  new  algorithm  to  limit  the  number  of 
interactions  between  N  vortex  elements  to  O(N),  instead  of  O(N^)  in  direct 
interactions,  using  multipole  expansion  of  the  contribution  of  groups  of 
elements. 

In  Appendix  I,  the  scheme  is  described  in  detail,  and  results  for  the 
evolution  of  a  temporal  shear  layer  are  analyzed.  A  temporal  shear  layer 
model  allows  one  to  limit  the  computations  to  a  fixed  number  of  large  eddies 
as  they  evolve  from  a  perturbation  to  a  coherent  vortex  structure.  The 
accuracy  of  the  computations  reveals  the  detail  of  the  inner  structure  of 
the  large  eddies  and  the  development  of  secondary  instabilities  which  force 
the  core  into  several  rotations,  as  well  as  the  severe  strain  which  material 
line  are  exposed  to  within  the  eddy  core.  The  application  of  the  scheme  to 
a  spatially  developing  shear  layer  is  presented  in  Appendix  II.  The 
improvement  of  the  accuracy  over  schemes  which  utilize  a  fixed  number  of 
elements  can  be  realized  by  comparing  these  results  with  results  shown  in 
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Appendix  III.  In  the  vortex  element  method,  the  number  of  elements  grows  as 
the  strain  field  develops,  letting  one  capture  the  areas  where  the  strain 
field,  and  thus  the  large  scalar  gradients  develop. 

Extension  of  the  vortex  element  method  to  three  dimensional  flow 
simulation  has  been  accomplished  by  utilizing  spherical  vortex  elements  that 
discretize  a  three  dimensional  vorticity  field  and  undergo  stretching  along 
the  direction  of  the  vorticity  vector,  as  described  in  Appendix  IV.  The 
application  of  the  scheme  to  study  the  evolution  of  azimuthal  instabilities 
on  vortex  rings  has  shown  that  these  structures  are  unstable  to  a  particular 
wave  number,  causing  the  ring  to  develop  into  a  star-like  structure  with 
lobes  of  vorticity  extending  in  the  radial  direction  of  the  ring.  Results 
also  show  that  these  azimuthal  instabilities  can  excite  higher  frequency 
modes  by  vortex  stretching,  generating  a  turbulent  cascade  of  the  energy 
into  higher  wave  numbers.  Results  of  the  simulation  compared  favorably  with 
experimental  data.  The  computations  are  currently  being  extended  to 
simulate  the  evolution  of  streamwise  vorticity  in  a  planar  shear  layer  and  a 
turbulent  jet. 

THE  TRANSPORT  ELEMENT  METHOD 

To  achieve  an  efficient,  self-adaptive  Lagrangian  algorithm  for  the 
solution  of  the  energy  and  species  conservation  equations,  scalar  gradients 
are  discretized  using  core  functions  similar  to  those  used  in  the  vortex 
element  method.  However,  contrary  to  the  scalar  concentration,  gradients 
are  not  conserved  along  a  particle  path  since  stretching  and  tilting 
material  layers  enhance  the  gradients  and  change  their  direction.  This 
effect  is  implemented  by  changing  the  strength  of  the  transport  elements 
according  the  variations  of  a  small  material  line  that  coincides  with  the 
center  of  the  transport  element.  Adding  a  chemical  source  term  requires 


changing  the  local  gradients  with  time.  This  amounts  to  transporting 
several  scalar  gradients  with  each  element  and  integrating  the  source  term 
within  each  element  to  compute  its  instantaneous  strength  according  to  a 
given  chemical  kinetics  scheme.  At  low  Mach  number,  volumetric  expansion 
due  to  heat  release  produces  an  irrotational  velocity  field,  and  generates 
non-baroclinic  vorticity  due  to  the  interaction  between  the  density  and 
pressure  gradients. 

In  Appendices  I  and  II,  the  mathematical  formulations  of  the  transport 
element  method  are  described  for  a  non  reacting  flow  and  a  for  reacting 
flow,  respectively.  In  Appendix  I,  the  method  has  been  applied  to  compute 
the  temperature  profiles  in  an  initially-thermally  stratified  temporal 
mixing  layer,  showing  how  entrainment  leads  to  intermittency  within  the  eddy 
core  and  to  mixing  enhancement  by  generating  large  gradients  as  the  material 
layers  stretch.  Statistics  of  mixing  of  a  passive  scalar  in  a  spatial  shear 
layer  have  been  compared  with  experimental  results  in  Appendix  III,  where 
only  a  simplified  version  of  the  scheme,  the  scalar  element  method,  was 
used.  The  comparison  is  favorable,  considering  the  fact  that  no  turbulence 
modeling  was  implemented  to  obtain  these  results.  The  accuracy  of  the 
computations  falls  off  around  the  boundaries  of  the  layer  due  to  the  small 
number  of  scalar  elements  which  were  used  in  this  simulation.  The  transport 
element  method  avoids  this  problem  by  transporting  the  gradients,  instead  of 
the  scalar,  and  is  expected  to  yield  better  predictions  for  the  mixing 
statistics.  This  is  currently  being  tested. 

Results  for  the  development  of  an  eddy  in  a  reacting  shear  layer  are 
presented  in  Appendix  II.  Initially,  reactants  and  products  are  on  the  top 
and  bottom  sides,  respectively.  As  the  eddy  grows  by  entraining  more 
reactants,  the  flame  is  stretched  and  wrinkled,  leading  to  a  rise  in  the 
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rate  of  product  formation  over  that  of  the  corresponding  laminar  flame. 
The  ratio  of  the  total  amount  of  products  formed  in  the  two  cases  scales 
with  the  square  of  the  flame  length  in  the  turbulent  case. 

To  validate  these  results,  computations  will  be  performed  for  a 
spatially  developing,  reacting  mixing  layer,  and  statistical  information 
derived  from  the  simulations  will  be  compared  with  experimental 
measurements . 

TURBULENCE-COMBUSTION  INTERACTIONS 

Slowing  down  the  rate  of  chemical  reaction  in  the  reacting  mixing  layer 
leads  to  local  and  temporary  extinction,  as  shown  in  Appendix  II.  As  the 
Damkohler  number  of  the  reacting  mixture  is  lowered,  the  rate  of  product 
formation  is  decreased.  Moreover,  the  reaction  is  observed  to  cease 
completely  for  short  periods  of  time.  The  lower  the  Damkohler  number,  the 
earlier  the  reaction  is  temporarily  extinct.  Local  extinction  is  observed 
around  areas  of  largest  strain  field.  The  extinction  occurs  temporarily 
since  on  the  other  side  of  the  layer,  products  at  high  temperature  heat  up 
the  reactants  and  resume  the  reaction  after  a  short  pause.  To  explain  why 
extinction  is  correlated  with  the  strain  field,  we  inspected  plots  of  the 
temperature,  strain  rate  and  expansion  rate  along  one  of  the  layers  within 
the  flame  zone.  It  was  found  that  a  negative  correlation  existed  between 
the  strain  rate  and  expansion  rate,  and  between  the  temperature  and  strain 
rate.  Thus,  it  was  concluded  that  as  the  strain  lowered  the  temperature  by 
enhancing  the  diffusion  flux  via  strong  gradients,  it  led  to  flame 
extinction.  The  temperature  drop  was  due  to  the  fact  that  chemistry  was 
slow  so  that  it  could  not  make  up  for  the  increase  in  the  diffusion  fluxes 
with  stretch. 
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Other  parameters,  e.g.,  the  Peclet  number,  Lewis  number  and  the 
activation  energy,  can  affect  the  interaction  between  turbulence  and 
combustion.  Work  is  underway  to  investigate  their  influence  on  the 
mechanism  described  above. 
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ABSTRACT 


In  computing  the  development  of  an  unstable  inviscid  shear  layer,  it  is 
found  that  using  a  fixed  number  of  vortex  elements  can  lead  to  large  errors 
due  to  the  strong  strain  field  which  develops  and  acts  to  distort  the 
original  vorticity  contours.  It  is  suggested  that  the  vorticity  should  be 
redistributed  among  elements  which  are  arranged  in  the  local  principal 
direction  of  strain  in  order  to  capture  this  distortion  accurately.  Mixing 
within  an  initially  stratified  layer,  which  results  from  the  combined  action 
of  convection  and  diffusion,  is  computed  using  a  similar  scheme  to  integrate 
the  energy  equation.  Calculations  illustrate  the  evolution  of  the 
temperature  profile  during  the  growth  of  the  instability. 


I  INTRODUCTION 


1.1.  BACKC3WUND 

Numerical  simulation  of  inviscid  two-dimensional  incompressible  flow 
using  vortex  discretization  of  the  Euler  equations  has  been  discussed 
extensively  in  recent  literature  (Leonard  [1],  Beale  and  Majda  [2],  Hald  [3] 
and  Ghoniem  and  Ng  [4]).  The  method  is  based  on  distributing  the  vorticity 
field  among  elements  which  carry  radially-symmetric ,  compact  supports  of 
vorticity  (Chorin  (5]).  By  choosing  the  extent  of  the  support,  or  the  core 
radius  of  each  element  to  be  larger  them  the  distance  of  separation  between 
neighboring  elements,  the  fields  of  individual  elements  overlap  and  high 
order  discretization  of  the  vorticity  field  can  be  achieved.  Vortex 
elements  move  with  the  local  flow  velocity  evaluated  at  their  geometrical 
centers,  which  is  computed  as  the  summation  over  the  contributions  of  all 
elements  that  exist  in  the  field.  The  motion  of  a  vortex  element  does  not 
change  its  circulation  and,  in  most  applications,  vortex  elements  possess 
invariable  core  shape  and  size. 

The  attraction  of  these  Lagrangian,  grid-free  methods  is  that,  by 
construction,  computational  vortex  elements  are  expected  to  be,  at  all 
times,  concentrated  around  zones  of  high  velocity  gradients.  When  properly 
exploited,  this  property  endows  the  scheme  with  the  resolution  necessary  to 
study  interesting  phenomena  that  arise  when  molecular  diffusion  is  small 
relative  to  convective  transport.  For  instance,  at  high  Reynolds  numbers, 
vorticity  exists  on  small  patches  of  the  fluid  and  it  suffices  to  distribute 
computational  elements  within  these  patches  and  hence  avoid  wasting  labor  on 
zones  of  very  small  vorticity.  That  the  elements  move  to  capture  large 
velocity  gradients  is  particularly  important  in  unsteady  and  nonlinearly 
unstable  flows  where  the  evolution  of  the  instability  causes  a  substantial 
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distortion  of  the  vorticity  distribution.  Moreover,  using  a  Lagrangian 
formulation  of  the  equations  of  motion  avoids  the  convective  non-linearity 
and  enables  the  construction  of  computational  schemes  which  are  explicit  in 
time.  The  employment  of  moving  Lagrangian  grids  (Fritts  and  Boris  [6]),  or 
grid-free  schemes  such  as  contour  dynamics  (Zabuski  et  al.  (7)),  are  other 
successful  ways  of  accomplishing  the  same  goal. 

1.2.  BRIEF  REVIEW 

Analysis  of  the  convergence  of  inviscid  vortex  methods  shows  that  three 
factors  govern  their  accuracy:  (1)  the  scheme  of  discretization  of  the 
initial  vorticity;  (2)  the  form  of  the  core  function;  and  (3)  the  ratio  of 
the  core  radius  to  the  separation  between  vortex  elements  (Chorin  et  al. 

[8],  Del-Prete  and  Hald  [9],  Hald  [3,10],  and  Beale  and  Majda  [2,11,12].) 
Results  of  these  analyses  have  been  supported  by  numerical  tests  (Nakamura 
et  al.  [13],  Roberts  [14],  and  Perlman  [15]).  In  the  following,  all  three 
factors  are  briefly  discussed. 

To  initialize  the  strength  of  vortex  elements,  Del-Prete  and  Hald  [9] 
used  the  average  vorticity  within  an  area  element  around  the  center  of  the 
element,  while  Beale  and  Majda  [2]  suggested  using  the  vorticity  at  the 
center  of  the  element.  Nakamura  et  al.  [13]  minimized  the  global  error 
between  the  continuous  and  the  discrete  vorticity  distribution  to  evaluate 
the  latter.  Anderson  and  Greengard  [16]  proposed  the  use  of  a  non  uniform 
mesh  to  discretize  the  vorticity  field.  Using  the  proceedure  in  [2]  or  [9], 
one  should  expect  almost  a  second-order  accuracy  for  short  time  if  the  core 
function  is  chosen  to  be  a  second  order  Gaussian.  A  fourth  order  Gaussian 
was  shown  to  improve  the  accuracy.  In  both  cases,  a  critical  parameter  is 
the  ratio  of  the  core  radius  to  the  distance  of  separation  between  the 
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centers  of  the  elements,  which  must  be  chosen  larger  than  unity  to  preserve 
the  accuracy  for  long  time. 

As  the  elements  move,  their  separation  exceeds  their  initial  value  if 
strong  strain  field  arises.  This,  in  effect,  decreases  the  critical  ratio 
of  core/separation,  leading  to  a  deterioration  of  the  accuracy.  The  fact 
that  large  strains  cause  deterioration  in  the  accuracy  of  vortex  methods  has 
been  observed  explicitly  in  analysis,  e.g.,  Leonard  [1].  Thus,  for  most 
inviscid  vortex  methods,  which  are  based  on  using  a  fixed  number  of  vortex 
elements  with  invariant  cores,  the  evolution  of  large  local  strains  can  lead 
to  large  errors.  For  example,  a  circular  patch  of  vorticity  may  deform  into 
an  elliptical  shape  with  its  major  axis  aligned  with  the  principal  direction 
of  strain.  If  a  small  fixed  number  of  computational  elements  is  used,  they 
may  not  be  able  to  accommodate  these  severe  changes.  Anderson  [17]  and 
Krasny  [18],  when  discretizing  non-smooth  vorticity,  employed  a  very  large 
core  radius  so  that  as  vortex  elements  moved  away  from  each  other  due  to 
stretch,  reasonable  overlap  could  still  be  maintained  to  satisfy  the 
requirements  for  accuracy.  One  may  also  be  forced  to  consider  schemes  of 
redistributing  the  vorticity  among  a  different  set  of  elements  under 
conditions  of  large  strain.  Similar  schemes  have  been  used  in  methods  of 
contour  dynamics  to  preserve  the  accuracy  of  the  integration  around  the 
vorticity  contours  (Zabuski  and  Overman  [19].)  Krasny  [20],  in  an 
independent  effort,  used  a  similar  procedure  in  simulating  the  evolution  of 
a  vortex  sheet  by  a  desingularized  Biot-Savart  integral. 

Extension  of  Lagrangian  element  methods  to  integrate  a  scalar 
conservation  equation  has  been  applied  to  several  problems  in  one  dimension 
(Chorin  [21],  Ghoniem  and  Oppeaheim  [22,23]  and  Ghoniem  and  Sherman  [24].) 
These  schemes  were  based  on  using  the  scalar  gradient,  in  analogy  to 
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vorticity,  in  the  transport  process.  Anderson  [16,25]  constructed  a  scheme 
to  solve  for  a  two  dimensional  thermal  in  the  inviscid  Boussinesq 
approximation  by  discretizing  the  density  equation  in  its  vortex  form.  This 
was  done  by  casting  the  equation  in  gradient  form  and  discretizing  the 
density  gradients  among  elements  that  could  be  transported.  This  scheme, 
while  preserving  the  advantages  of  the  vortex  method,  suffers  from  a  major 
problem:  A  large  strain  field,  while  it  may  lead  to  the  generation  of  large 
gradients,  depletes  the  area  of  computational  elements  which  are  used  to 
transport  these  gradients. 

1.3.  ORGANIZATION 

In  this  paper,  we  apply  the  inviscid  vortex  methods  to  the  problem  of  a 
temp  ral  shear  layer  at  high  Reynolds  number.  This  problem  is  characterized 
by  a  well-defined  smooth  vorticity  field  at  time  zero,  and  has  well- 
documented  stability  properties.  At  later  times,  the  shear  layer  develops 
into  a  complicated  structure  which  resembles  a  turbulent  eddy,  and  a  very 
strong  strain  field  is  generated  around  this  eddy.  We  use  the  analytical 
solution  of  a  temporal  shear  layer  to  measure  the  accuarcy  of  the  results  at 
the  initial  stages  of  development,  and  test  the  schemes  for  initializing  the 
vortex  elements.  At  longer  times,  we  observe  the  effect  of  the  strain  field 
on  the  accuracy  of  the  computations  and  suggest  ways  to  cope  with  it.  We 
then  proceed  to  compute  the  temperature  field  as  fluids  with  different 
temperatures  are  entrained,  stretched  and  mixed. 

In  Section  II,  the  formulation  of  the  vortex  method  is  described,  and 
is  extended  to  solve  for  a  flow  with  a  strong  strain  field.  The  scheme  is 
applied  to  compute  the  evolution  of  a  vorticity  layer  subject  to  periodic 
boundary  conditions.  The  growth  of  the  instability  and  its  effects  on  the 
flow  field  are  investigated,  in  Section  III,  the  concepts  of  the  vortex 


method  are  generalized  to  solve  the  energy  equation  and  to  obtain  the 
temperature  profile  across  the  shear  layer  during  its  development.  The 
paper  ends  with  conclusions  in  Section  IV. 


II  INVISCID  INCOMPRESSIBLE  FLOW 
II. 1.  THE  VORTEX  METHOD 

For  am  inviscid  incompressible  flow,  the  vortex  transport  equation  is: 

+  u-Vw  -  0  (1) 

4  ♦  ■  *  «  (2) 

where  u  -  (u,v)  is  the  velocity,  w  -  V  x  u  is  the  vorticity,  x  -  (x,y)  are 
the  strean&rise  and  cross  stream  directions,  respectively,  t  is  time,  V  - 
(3/3x,3/3y)  and  A  -  V.V.  Variables  are  normalized  with  respect  to  the 
appropriate  combination  of  a  characteristic  velocity  and  length  scale.  *  is 
the  stream  function  defined  so  that  u  -  3^/3y  and  v  -  -  ty/dx.  The  solution 
of  Eq.  (1)  can  be  written  as: 

oo<X(X,t),t)  -  «(X,0)  (3) 

while  x  is  governed  by: 

^  -  u(X(X,t),t)  (4) 

where  X(X,0)  -  X.  In  the  vortex  method,  the  vorticity  field  w(X,0)  is 
discretized  between  elements  centered  at  X^  i-l,..,N,  so  that: 


<*>(x,0)  -  Z  T.  f.{x-X. )  (5) 

i-1  1  6  l 

2 

where  h  is  the  circulation  of  an  element  of  strength  and  f^  is  the 

core  function.  fj(x)  -  1/62  f(r/S),  where  r2-x2+y2,  and  J  f&  dx  -  1.  S  is 
the  core  radius,  and  f j  is  a  fast  decaying  function  so  that  most  of  the 
vorticity  is  concentrated  within  r  <  5.  To  approximate  the  initial 
vorticity  distribution  accurately,  5  should  be  greater  than  h,  where  h  is 
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the  initial  separation  between  vortex  centers.  The  core  function  f  plays  a 
similar  role  as  interpolating  polynomials  in  finite-difference  schemes  and 
base  functions  in  finite-element  formulations.  By  requiring  f  to  be 
radially  symmetric,  the  approximation  in  Eq.  (5)  is  at  least  second  order. 

Using  Eq.  (3)  and  the  incompressibility  condition,  the  vorticity 
distribution  at  any  time  is  given  by 


N 

w(x,t)  -  E  T.  f c(x-X:  )  (6) 

i-1  1  S  ^ 

where  dxi/dt-u(xi>t)  and  x^X^OJ-OL  . 

The  stream  function  of  a  single  vortex  element  is  obtained  by 
integrating  Eq.  (2).  Using  polar  coordinates,  for  a  vortex  element  placed 
at  x-0,  3*5/3r  -  -x(r/4)/r,  where  tc(r)  -  QJr  r'  f(r')  dr'.  Moreover,  u0-  - 
3tj/3r.  The  velocity  field  induced  by  a  distribution  of  vortex  elements  ,  of 
shape  f&  and  strength  I\  located  at  Xi(Xi,t)  is: 

J  ri  I's(x-Xi>  (7) 


where 


(8) 


Vortex  elements  move  without  changing  their  circulation  (strength)  or  core 
shape,  at  a  velocity  computed  from  Eq.  (7). 

In  the  calculations,  we  used  mostly  a  second  order  Gaussian  core: 


When  applying  the  vortex  scheme  to  a  flow  field  with  boundary 
conditions  other  than  u(«,t)«0,  a  potential  flow  is  added  to  satisfy  these 
conditions.  In  this  work,  we  perform  computations  for  a  periodic  shear 
layer.  The  velocity  field  induced  by  vorticity  outside  the  computational 
domain  0  <  x  <  X,  where  X  is  the  longest  wavelength  of  the  perturbation, 
must  be  added  to  u^.  The  total  velocity  is: 


“  Fi 

U  «  l  y- 

i-1  2n 


i  f  <r -  -(xi3xl>  exp  <- 

j-0  ((x+jXr+y^) 


((x+jX)2+y2) 


.  n  (~sinh(2ny/X) ,  sin(2nx/X))  , 

\  ( cosh( 2ny/X)-  cos( 2nx/X) ) 

where  N  is  the  total  number  of  vortex  elements  in  the  computational  domain  0 
<  x  <  X.  Note  that  since  S  <<  X,  the  effect  of  the  core  was  included  only 
for  the  nearest  sister  vortices. 

The  initial  vorticity  distribution  across  the  shear  layer  can  be  well 
represented  by  a  Gaussian  curve  (which  should  not  be  confused  with  the 
Gaussian  core  of  individual  vortex  elements)  with  a  spread  2 a: 


2(X)  -  —  exp(-  Y2 /a2 )  (11a) 

/A.<r 

where  AU  is  the  velocity  difference  across  the  layer  and  a  is  the  standard 
deviation  of  the  Gaussian.  The  corresponding  velocity  distribution  is: 


U(X)  -  ip  erf(Y/«)  (lib) 

x  2 

where  erf(x)  -  2//A.  gj  exp  (-  r  )  dr  is  the  error  function.  We  take  AU  and 
a  as  the  characteristic  velocity  and  length  scales  of  the  problem, 
respectively. 
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As  it  was  pointed  out  in  the  Introduction,  using  either  a  pointwise 
discretization,  or  an  area  average  value  2(X)  dx,  where 

Xj,  i-1,2,  ..  ,N  are  the  centers  of  a  square  mesh  of  side  h,  to  discretize 
the  vorticity  of  the  shear  layer  among  vortex  elements  produced  a  large 
error  in  the  initial  growth  of  the  perturbation.  Instead,  the  following 
scheme  was  used: 

N  2 

fi(X.  )  -  £  w.  hZ  f.(X.-X.)  (12) 

i  j.l  1  6  1  ^ 

for  i-l,2,...,N.  The  error  associated  with  this  distribution  was  used  as  a 
measure  of  the  accuracy  of  the  initial  discretization.  In  all  cases,  the 
error  )ew)  -  J  |Q(X)-<»>(X,0)  |  dx  <  10~^.  The  error  ew  increased  rapidly  as 
S/h  was  decreased  below  one,  which  is  consistent  with  the  result  of  the 
convergence  theory  which  shows  that  the  overlap  between  neighboring  elements 
is  necessary  for  accurate  discretization  of  vorticity.  For  S/h  >  1,  the 
error  was  less  sensitive  to  its  exact  value,  uptil  S/h  -1.5.  In  the 
following  calculations,  we  used  S/h  -  1.1  -  1.4. 

To  measure  the  effect  of  the  accuracy  of  the  initial  discretization  of 
vorticity  among  vortex  elements  on  the  flow  field  for  short  times,  we  will 
use  the  rate  of  growth  of  the  perturbation.  The  growth  of  the  initial 
perturbation  can  be  characterized  by  an  integral  parameter  I  as: 

I  -  0;X_J-  |u(x,t)  -  U(x)  |  dx  (13) 

which  is  used  in  the  linear  theory  analysis  of  the  perturbation. 

At  t  -  0,  the  layer  was  perturbed  by  a  sinewave  with  amplitude  e,  taken 
as  0.001  X,  0.01  X,  and  0.1  X.  In  Fig.  1,  we  compare  the  growth  of  the 
perturbation  with  the  prediction  of  the  linear  theory  of  stability  (Michalke 
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[26].)  to  assess  the  accuracy  of  the  vortex  method  for  short  times.  For 
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most  of  the  computations,  X  -  X  -  13.2  a,  which  corresponds  to  the  wave 
with  the  maximum  growth  rate.  Eq.  (4)  is  integrated  using  a  second  order 
Heun's  method  with  at  -  0.1,  and  h-X/44-0.3,  and  N(0)  -  572  vortex 
elements.  The  figure  indicates  that  for  e  -  0.01  and  0.001,  the  layer 
behaves  linearly  and  the  computed  growth  rate  i  -  dlnl/dt  -  0.215  agrees 
well  with  the  results  of  the  linear  theory,  i  -  0.22.  The  latter  was 
computed  as  the  eigenvalue  of  the  linearized  Euler  equations  (Betchov  and 
Criminale  [27].)  Using  h  -  X  /  24,  i.e.  N(0)  -  168,  and  a  second  order  time 
integration  scheme,  i  -  0.23.  For  N(0)  -  572  and  a  first  order  time 
integration,  i  -  0.24.  Within  this  linear  stage  of  development,  the  maximum 
distance  between  neighboring  element  in  the  direction  of  maximum  strain  Ax  < 
1.5  h,  i.e  the  flow  is  developing  mild  stretch.  For  e  -  0.1,  the 
perturbation  leads  directly  to  the  nonlinear  range. 

In  Figs.  2,  3  and  4,  the  vortex  elements  and  their  velocity  vectors  are 
plotted  for  e  -  0.001  X,  0.01  X,  and  0.1  X,  respectively.  In  the  first  two 
cases,  the  end  of  the  linear  range  corresponds  to  the  beginning  of  the 
rollup  of  the  interface,  defined  here  as  the  line  which  coincides  with  y  -  0 
at  t  -  0,  and  the  formation  of  a  spiral  center  at  the  midpoint  of  the 
wavelength.  Concomitantly,  the  interface  starts  to  stretch  near  the 
boundaries  of  the  domain  and  two  saddle  points  are  established  at  the 
beginning  and  end  of  the  wavelength,  x  ■  0,  and  X.  Beyond  the  linear  range, 
the  perturbation  continues  to  grow  with  more  layers  rolling  around  the 
spiral  center  and  stretching  near  the  saddles.  Within  this  nonlinear  range 
of  development,  special  care  must  be  exercised  or  the  numerical  accuracy 
deteriorates  quickly,  as  exhibited  by  the  evolution  of  irregular  motion  near 
the  saddles  and  the  loss  of  organization  of  the  evolving  structure. 


II. 2.  EFFECT  OF  STRETCH 

The  loss  of  organization,  which  is  associated  with  the  development  of 
strong  stretch,  illustrates  one  of  the  fundamental  problems  of  the  vortex 
method.  Vortex  elements,  which  start  as  cores  with  radial  symmetry,  may  not 
properly  represent  the  vorticity  field  after  it  has  developed  strong  local 
strains.  As  the  effective  distance,  Ax,  between  neighboring  elements 
increases,  the  ratio  6/Ax  (equivalent  to  S/h)  reaches  levels  where  the 
vorticity  discretization  becomes  inaccurate.  One  obvious  remedy  is  to 
restart  the  calculations  with  smaller  values  of  h  to  allow  a  larger  number 
of  weaker  elements  to  represent  the  strong  distortion.  However,  that  only 
delays  the  onset  of  the  crisis  at  the  expense  of  using  more  elements  at  the 
initial  stages  when  they  are  not  needed.  Several  remedies  may  be  suggested: 
(1)  utilizing  deformable  cores;  (2)  employing  large  cores;  or  (3)  using  more 
elements  as  the  distance  between  the  original  elements  increases. 

The  first  scheme,  utilizing  deformable  cores,  depends  on  assuming  that 
the  core  structure  will  become  elliptical  as  stretch  develops,  with  the 
major  axis  of  each  element  aligned  with  the  local  principal  direction  of 
strain.  The  vorticity  distribution  within  the  core  must  also  adapt  to  the 
geometrical  boundaries  of  the  cores.  If  elements  with  constant  vorticity 
within  the  cores  and  zero  outside,  i.e.  Rankine  vortex  elements,  are  used 
then  these  elements  will  become  Kirkchoff  vortices  which  have  analytical 
expressions  for  the  induced  velocity  field.  However,  there  is  an  obvious 
limitation  on  maintaining  one  ellipse  as  a  single  element  if  the  ratio 
between  its  axes  exceeds  a  reasonable  value.  Thus,  this  scheme  is 
discarded. 


The  second  scheme,  in  which  one  uses  large  cores,  did  not  yield 
accurate  predictions  for  the  growth  rate  within  the  linear  range,  in 
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accordance  with  the  results  of  the  convergence  theory.  Moreover,  it  will 


fail  at  the  point  where  6/AX  <  1  due  to  stretch.  It  does,  however,  delay 


the  deterioration  of  accuracy  since  it  maintains  a  reasonable  overlap 


between  neighboring  elements  for  longer  times. 


The  third  option,  redistributing  the  vorticity  field  among  an 


increasing  number  of  elements  arranged  along  the  direction  of  principal 


direction  of  strain,  is  employed  here.  One  monitors  the  distance  between 


neighboring  elements  in  the  direction  of  maximum  positive  stretch  AX.  If  Ax 


>  p  h,  where  1  <  0  <  2,  an  extra  element  is  placed  halfway  between  the 


original  elements  and  the  vorticity  is  redistributed  to  compute  the  share  of 


the  new  element.  Ideally,  this  redistribution  should  not  perturb  the 


existing  vorticity  field,  that  is 


N  N+n  ~ 

rt  fjU-V  -  rt  fjU-V 


where  n  is  the  number  of  new  particles,  and  a  ~  indicates  the  new  value  of 


the  strength  and  location  of  the  vortex  elements.  Unfortunately,  this  is  a 


large  dense  system  of  linear  equations  to  be  solved  every  time  step. 


Therefore,  its  benefit  does  not  warrant  the  added  cost. 


A  more  economical  scheme  is  based  on  equally  interpolating  the  strength 


of  the  two  original  elements  among  the  three  elements,  i.e.  assuming  uniform 


stretch  between  the  two  original  elements.  This  amounts  to  splitting  the 


original  vortex  dumbbell  formed  of  two  vortex  discs  into  three  discs  when 


the  distance  between  the  centers  of  the  two  discs  exceeds  a  threshold,  as 


shown  in  Fig.  5.  To  minimize  the  interpolation  errors,  the  maximum 
interdistance  between  neighboring  elements  is  taken  as  1.5  h.  This  will 


also  keep  the  ratio  5/Ax  within  reasonable  limits. 
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To  illustrate  the  degree  of  stretch  experienced  by  this  flow,  we  plot 
the  growth  of  the  length  of  the  interface,  and  the  total  number  of  vortex 
elements,  N(t),  used  to  capture  this  stretch  for  three  perturbations  c  - 
0.001X,  0.01X,  and  0.1  X  in  Figs.  6,  7,  respectively.  Within  the  linear 
range  the  layer  is  subjected  to  mild  stretch  and  N  remains  almost  constant. 
Beyond  that,  the  length  of  the  line  grows  linearly  and  N  multiplies 
accordingly.  From  the  plots  of  the  location  of  vortex  elements,  we  noticed 
that  most  of  the  stretch  is  concentrated  around  the  spiral  center  and  the 
saddles  at  the  boundaries  of  the  domain. 

II. 3.  SHEAR  LAYER  DYNAMICS 

Figures  1,  2,  3  and  4  reveal  that  the  growth  of  the  perturbation  and 
the  development  of  the  eddy  structure  can  be  divided  into  four  stages:  (1) 
linear  growth;  (2)  rise  to  a  maximum  amplitude;  (3)  decay  to  a  constant 
amplitude;  and,  (4)  very  slow  decrease  of  amplitude.  The  first  stage  has 
been  discussed.  The  strongest  stretch  and  fastest  multiplication  of  the 
vortex  elements  occur  during  the  second  stage  where  an  eddy  is  forming  in 
the  middle  of  the  wavelength  and  two  braids  are  evolving  between  each  two 
neighboring  eddies.  During  this  stage,  the  core  maintains  almost  a  circular 
configuration  and  the  stretch  is  concentrated  within  the  braids. 

In  the  third  stage,  the  eddy  deforms  into  an  elliptical  structure, 
while  the  size  of  the  perturbation  decreases  from  its  maximum  value.  This 
is  accompanied  by  more  stretch  along  the  braids  and  within  the  core,  and  a 
slovxiown  of  the  eddy  rotation.  By  the  end  of  this  stage,  the  thickness  of 
the  braids  at  the  saddle  points  has  become  extremely  small.  At  the  final 
stage,  the  envelope  of  the  core  reaches  a  dynamic  equilibrium,  i.e.,  it  does 
not  rotate  any  more,  while  its  boundaries  keep  on  stretching  as  the  fluid 
within  the  eddy  starts  to  move  in  the  main  directions  of  the  streams. 
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Although  there  are  signs  of  that,  it  is  difficult  to  confirm  that  the  flow 
has  reached  a  steady  state. 

The  kinetic  energy  of  perturbation  u' .  u'/2,  where  u'-  u  -  U,  and  the 
total  kinetic  energy  in  the  flow  within  the  computational  domain,  \i.\i/2,  are 
plotted  in  Figs.  8a  and  8b.  The  first  quantity  rises  with  the  growth  of  the 
perturbation  and  the  formation  of  the  eddy,  then  falls  with  the  collapse  of 
the  eddy  and  the  return  of  the  fluid  to  the  main  streams  (Corcos  and  Sherman 
[28]).  The  total  kinetic  energy  is  conserved  since  the  flow  is  inviscid. 

Using  larger  values  for  h  while  keeping  S/h  the  same  caused  a  slight 
fattening  of  the  core  at  latter  times,  while  the  main  features  of  the  flow 
were  reproduced  almost  exactly.  A  similar  modification  of  the  structure  is 
observed  when  using  a  first  order  time  integration  scheme,  or  increasing  the 
time  step.  It  was  concluded  that  the  errors  introduced  by  using  a  small 
number  of  elements  or  a  low  order  time  integration  scheme  were  numerical- 
diffusion-like  errors.  We  also  found  that  the  dependence  on  the  value  of  h, 
or  the  initial  number  of  elements,  becomes  much  less  pronounced  when  the 


scheme  of  increasing  the  number  of  elements  with  stretch  is  employed. 


Figure  9  shows  a  qualitative  comparison  between  the  experimental  results  of 


Roberts  et  al.  [29]  and  the  computational  results.  Here  we  use  a  Galilean 


transformation  to  compare  the  experimental  results  of  the  spatially- 


developing  layer  and  the  computational  results  of  the  temporal  layer. 


The  physical  parameters  that  govern  the  flow  field  are  X  and  e. 


Results  for  the  rollup  of  a  layer  with  X  -  10.5  <  X  are  presented  in  Figs. 


10  and  11,  showing  the  growth  of  the  perturbation  and  the  vorticity  field. 


X  is  the  wavelength  of  the  most  unstable  perturbation.  The  computed  growth 


rate  I  -  0.214  while  the  analytical  value  is  0.208.  More  vorticity  remains 
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in  the  braids  between  the  eddies  which  are  not  strong  enough  to  accomplish 
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the  same  stretch  as  in  the  case  of  X  . 

Figures  12  and  13  show  results  for  X  -  2  X*,  with  e  -  0.01X  and  0.1X, 
respectively.  In  the  first  case,  i  -  0.18  while  the  analytical  result  is 
0.173.  The  core  is  smaller  and  weaker  than  for  the  case  of  X*  and  hence  the 
braids  are  thicker  and  maintain  more  of  the  original  vorticity.  At  later 
times,  a  small  scale  rollup  is  observed  near  the  boundary  of  the  domain  due 
to  the  instability  of  the  vorticity  layer  that  forms  the  braids.  This 
rollup  occurs  only  at  the  fourth  stage  of  development  when  the  midsection  of 
the  braids  becomes  almost  stationary,  i.e.  when  the  motion  produced  by  the 
braids  is  neutralized.  Comparing  Figs.  12  and  13,  we  see  that  contrary  to 
the  most  unstable  case,  the  effect  of  the  initial  perturbation  is  more 
pronounced  here  in  terms  of  the  size  and  shape  of  the  eddy  and  the  braids. 
Higher  amplitudes  of  perturbation  tend  to  form  a  larger  core  and  thinner 
braids.  The  ratio  between  the  major  and  minor  axes  of  the  elliptical  core 
increases  with  e  and  small  amplitude  waves  start  to  appear  on  the  braids. 

Figures  14  and  15  show  results  for  X  -  3  X*  with  amplitudes  c  «  0.01X, 
and  0.1X,  respectively.  The  effect  of  the  amplitude  is  emphasized  further 
since  at  larger  c,  the  core  splits  into  two  eddies.  This  bifurcation 
phenomenon  was  observed  before  by  Pozrikidis  and  Higdon  (30).  The  braid 
instability  is  manifested  here  by  the  long  waves  that  appear  at  the  later 
stages  of  development  of  the  layer. 

With  the  presence  of  two  perturbation  wavelengths,  a  new  process  is 

observed.  Figures  16  and  17  depict  results  for  a  layer  subject  to  two 

★  ★  ★ 

perturbations  superimposed  at  t  -  0,  at  X  and  2X  with  c  -  0.1X  for  both 
perturbations.  The  results  show  that  when  the  amplitude  of  the  two 
perturbations  are  equal,  pairing  starts  at  the  end  of  the  second  stage  and 
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before  any  substantial  elongation  of  the  eddies.  The  growth  of  the 
subharmonic  perturbation  closely  resembles  that  of  the  fundamental,  as  shown 
in  Fig.  17.  The  eddies  continue  to  deform  while  they  pair  until  the  "vortex 
fluid"  contained  within  each  structure  start  to  rotate  around  a  common 
center  and  their  original  boundaries  become  indistinguishable.  Similar 
qualitative  observations  were  shown  in  the  computations  of  Corcos  and 
Sherman  (28]  and  Riley  and  Metcalfe  (31). 


Ill  THE  TEMPERATURE  DISTRIBUTION 


III.l.  THE  TRANSPORT  ELEMENT  METHOD 

In  an  inviscid  incompressible  flow,  the  temperature  distribution 
evolves  according  to  the  following  form  of  the  conservation  of  energy: 

+  u* VI  -  0  (15) 

where  T  is  temperature.  This  is  equivalent  to  the  statement  that 
T(x(*»t),t)  -  T(  X,  0 ) ,  where  x(X,0)  -  X  and  dx/dt  -  u(X,t).  To  solve  this 
equation  using  a  Lagrangian  element  scheme,  we  start  by  introducing  the 
temperature  gradient  q  -  VT,  where  q  -  (p,q)  is  a  vector  proportional  but 
opposite  to  the  heat  flux  vector  -k  q,  k  being  the  thermal  conductivity. 

The  transport  equation  of  q  is  obtained  by  taking  the  gradient  of  Eq.  (15) 
and  rearranging: 

J9  +  u*7q  -  -  q.9u  -  q  x  «  (16) 

where  <•>  -  w  e„  and  e_  is  the  unit  vector  normal  to  the  (x,y)  plane.  Thus, 
along  a  particle  path  X(X,t),  the  temperature  gradient  changes  according  to 
the  local  strain  field  and  turns  with  the  local  rotation  of  the  fluid 
element.  Using  the  vortex  method  described  in  the  previous  section,  the 
velocity  gradient  may  be  computed  directly  from  the  vorticity  distribution 
as:  Vu  -  E  r.  VK.(x-X:  )  +  Vu,  where  u  is  the  irrotational  component  of  the 

1  0  1  p  p 

velocity. 

The  scheme  proceeds  in  the  same  way  as  the  vortex  algorithm.  The 
initial  temperature  gradient  is  discretized  among  a  number  of  elements 
located  at  the  centers  of  a  square  mesh  of  side  h  so  that: 
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q(x,0) 


N  2 

£  Qi  h  f.(x-X.) 
i-1  1  61 


(17) 


where  q(x,0)  is  the  initial  distribution  of  the  temperature  gradient,  q  - 
VT(X,0).  To  initialize  a  similar  procedure  to  that  used  in  computing 
the  strength  of  the  vortex  elements  is  employed  here,  i.e.  Eq.  (12)  with  q 
instead  of  w,  is  solved  to  find  qi(0).  To  update  qi ( t ) ,  Eq.  (16)  is  solved 
in  two  fractional  steps:  in  the  first  step,  the  elements  are  transported 
without  changing  their  strength  or  their  core  shape  or  size.  In  the  second 
step,  the  strength  of  the  elements  is  updated  according  to: 

dq. 

ar  -  -  qt  •  ^i  -  X  «i  us) 

Thus,  a  system  of  ordinary  differential  equations  must  be  integrated  to 
update  the  strength  of  the  gradient  elements  as  they  move  along  particle 
paths.  The  local  gradient  at  time  t  is  computed  from: 


q(x,t) 


t  qt(t)  h  f6(x-xi) 


(19) 


The  core  function  f  ^  may  be  different  for  the  vortex  elements  and  the 
gradient  transport  elements.  In  this  work,  we  use  the  same  form  for  both. 
The  temperature  can  be  calculated  by  direct  integration  of  the  gradient 
along  a  determined  path.  As  pointed  out  by  Anderson  [17],  a  convenient 
expression  can  be  obtained  by  expressing  the  temperature  as  a  Poisson 
integral  in  the  temperature  gradient,  T  -  J  73(x-x' )  .  VT(x-x')  dx' .  Using 
Eq.  (19)  for  q^,  we  get: 


N 

T(x,t)  -  E^  qt(t)  .  VG^x-jq)  (20) 
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while  the  relationship  between  f  and  k  is  as  before.  This  expression  is 
convenient  when  used  in  connection  with  the  vortex  method  since  all  the 
functions  involved  in  the  summation  must  be  computed  for  the  transport  of 
the  vortex  elements  and  with  simple  programming  tricks,  the  increase  in  cost 
can  be  minimized. 

Results  obtained  for  the  temperature  distribution  in  the  shear  layer  of 
Figs.  2,  3  and  4  are  shown  in  Figs.  18,  19  and  20,  respectively.  At  t  -  0, 
the  temperature  distribution  is  described  by  an  error  function,  with  T(x,0) 

-  0.5  (1  +  erf(Y)).  This  choice  is  motivated  by  the  fact  that  this  is  the 
fundamental  solution  of  the  diffusion  equation.  Therefore,  an  initial 
discontinuity  in  temperature  would  develop  into  an  error  function  before  the 
perturbation  grows  and  convection  effects  become  important.  The  layer  is 
first  perturbed  by  a  sinewave  by  displacing  the  elements  according  to 
Y-Y(X),  and  then  the  temperature  gradient  Q(X)  is  computed.  The  discrete 
values  qi(0)  are  obtained  as  follows:  Since  the  temperature  is  constant 
along  the  streamlines  after  the  perturbation  Y-Y(X),  then  T(X,0)  -  0(X)  - 
0.5(1  +  erf(Y  -  Y(X) ) .  From  this  equation  we  can  recover  the  initial 
distribution  of  P(X,0),  and  Q(X,0)  as  follows: 


P(X)  -  ||  «  -Gau(  Y-Y(X) )  gj- 
Q(X)  -  ||  -  Gau( Y-Y(X) ) 


(21) 


20 


where  Gau  is  a  Gaussian  similar  to  Eq.  (11a).  These  expressions  are  then 
used  in  Eq.  (  17)  to  compute  qi-  The  values  of  q^(0)  were  initialized  one 
column  at  a  time,  i.e.  for  fixed  values  of  X,  to  avoid  solving  N2 
simultaneous  equations,  and  instead  solve  Nx  of  simultaneous  equations. 
The  error  associated  with  this  approximation  was  very  small  since  the 
perturbation  was  kept  at  a  low  value. 

In  the  computations,  we  used  the  same  particles  to  transport  vortex 
elements  and  elements  of  the  temperature  gradient.  This  represents  a 
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substantial  saving  since  the  kernel  functions  appearing  in  the  expressions 
of  the  velocity,  velocity  gradients,  and  temperature  can  be  computed  all  at 
once  and  the  velocity  is  computed  only  for  one  set  of  elements. 

III. 2.  ENTRAINMENT  IN  A  SHEAR  LAYER 

To  quantify  the  overall  change  in  the  temperature  distribution,  we 
define  a  quantity  T,  similar  to  the  growth  I,  as; 

T  -  0;X_j"  |T(x,t)  -  e(x)|  (22) 

T  can  be  regarded  as  an  average  thermal  thickness  of  the  shear  layer. 

Within  the  linear  range,  the  temperature  distribution  remains  essentially 
the  same,  except  for  getting  shifted  up  or  down  depending  on  the  local  sign 
of  the  perturbation.  In  Fig.  21,  the  natural  logarithm  of  T(t)  is  shown  for 
three  values  of  the  initial  amplitude  of  the  perturbation.  The  accuracy  of 


the  calculation  of  the  temperature  profiles  depends  on  the  initialization  of 
the  vorticity  and  temperature  gradient,  and  on  the  value  of  6/h. 

During  the  second  stage,  and  with  the  rollup  of  the  interface  and  the 
establishment  of  a  spiral  center  at  the  midpoint  of  the  wavelength,  a 
complex  temperature  gradient  develops  as  a  result  of  the  motion  of  the  cold 
fluid  upwards  and  the  hot  fluid  downwards  around  the  spiral  center.  Within 
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this  stage,  if  the  number  of  transport  elements  remained  the  same,  i.e. 
stretch  was  not  accommodated  by  introducing  elements  where  the  local  strain 
is  large,  the  temperature  distribution  would  collapse  very  quickly.  In  the 
problem  of  periodic  shear  layer,  this  collapse  leads  to  values  of  T(x,+-,t) 

<  1  and  T(x,-«,t)  >  0.  The  reason  of  the  loss  of  accuracy  is  clear  from  Eq. 
(16).  When  the  elements  move  apart,  the  accuracy  of  computing  the  velocity 
gradient  Vu  deteriorates,  and  hence  the  new  values  of  q^  accumulate  large 
errors.  Thus,  while  the  calculation  of  the  velocity  field  at  the  early 
stages  of  strong  stretch  using  a  fixed  number  of  vortex  elements  may  be 
acceptable  for  a  short  period  of  time,  the  calculations  of  the  velocity 
gradient  and  the  evolution  of  a  passive  scalar  will  show  unacceptable 
errors. 

To  continue  beyond  the  linear  stage,  the  distance  between  neighboring 
elements  in  the  principal  direction  of  strain,  Ax,  must  to  be  monitored.  If 
Ax  >  (3h,  where  0  >  1,  one  extra  element  is  added  between  the  two  original 
elements  and  the  total  value  of  q^  is  redistributed  equally  between  the 
three  elements.  In  the  calculations,  we  used  0  -  1.5.  Numerical 
convergence,  in  which  one  systematically  refines  the  numerical  parameters 
until  no  more  changes  are  observed,  was  used  to  obtain  these  results. 

The  effect  of  the  shear  layer  rollup  on  the  temperature  distribution  is 
seen  in  Figs.  18  end  19.  Immediately  after  the  interface  reaches  a  vertical 
position,  an  S-shape  starts  to  form  indicating  that  cold  fluid  has  been 
transported  from  the  lower  stream  into  the  upper  stream  and  vice  versa. 

This  phenomenon,  known  physically  as  engulfment  or  entrainment,  relies 
solely  on  convective  transport  and  is  observed  when  molecular  diffusion, 
which  acts  to  dissipate  the  sharp  gradients,  is  small.  Fast  entrainment 
with  small  diffusion  leads  to  "unmixedness"  of  hot  and  cold  fluids  within 


the  eddy  core.  With  more  fluid  being  transported  to  the  opposite  stream  the 
S- shape  grows,  reaching  a  maximum  amplitude  when  the  interface  becomes 
horizontal.  At  this  moment,  fluid  with  the  maximum  and  the  minimum 
temperature  has  been  entrained  into  the  core,  i.e.  entrainment  has  reached 
all  the  way  to  the  free  streams  to  bring  fluid  into  the  core  of  the  eddy. 
This  is  the  stage  of  maximum  entrainment  when  the  core  size  reaches  its 
largest  size  and  cannot  accommodate  any  more  fluid.  In  the  case  of  c  - 
0.1X,  it  corresponds  to  t  -  8.0,  which  is  at  the  end  of  the  second  stage  of 
development.  To  make  the  correspondence  between  the  temperature  profiles 
and  the  evolution  of  the  interface  of  the  layer  clear,  we  plot  the  latter  in 
Fig.  22,  showing  the  actual  vortex  elements  that  were  used  in  the 
computations  of  the  interface.  At  this  time,  the  interface  has  rotated  180° 
around  the  spiral  center.  This  is  the  first  step  in  the  process  of 
homogenization  of  the  core. 

As  the  core  rotates  further  into  the  third  stage,  the  inner  part  of  the 
interface  develops  a  secondary  instability  that  rolls  up  in  a  very  similar 
manner  to  the  primary  instability.  This  secondary  instability  is  in  phase 
with  the  primary  instability  and  can  be  envisioned  by  zooming  in  on  the 
intersection  between  the  interface  and  the  horizontal  centerline  of  the 
layer.  Due  to  the  elongation  of  the  outside  envelope  of  the  core,  the 
wavelength  of  the  secondary  instability  grows  with  time,  as  seen  from  Fig. 
22.  However,  the  amount  of  fluid  within  the  elliptical  envelope  remains 
constant,  or  decreases  slowly  as  seen  from  Fig.  21  for  the  temperature 
thickness  of  the  layer.  The  growth  of  the  secondary  instability  provides  a 
mechanism  of  internal  entrainment  within  the  core.  During  the  growth  of  the 
secondary  instability,  an  inverted  S-shape,  or  a  Z-shape,  forms  in  the 
middle  of  the  temperature  profile,  Figs.  18-20.  The  entrainment  associated 


with  this  instability  turns  the  fluid  in  a  clockwise  fashion,  making  the 
inside  of  the  core  more  uniform.  This  is  seen  from  the  decay  of  the  peaks 
in  the  temperature  profile  as  this  Z-shape  grows. 

With  another  180°  turn  of  the  interface  at  the  spiral  center,  a  smaller 
S-shape  forms  in  the  middle  of  the  profile  due  to  '.he  onset  of  an  even 
shorter  wavelength  instability  that  is  in-phase  with  the  primary 
instability.  While  the  existence  of  the  secondary  instability  was  not 
observed  before  in  numerical  simulations,  its  presence  is  clearly  seen  in 
the  experimental  results  in  Fig.  9. 

The  onset  and  subsequent  growth  of  successively  shorter  wavelength 
instabilities  continues,  leading  to  a  more  uniform  temperature  distribution 
within  the  eddy  core.  An  asymptotic  limit  to  this  process  can  be  foreseen: 
it  is  the  formation  of  a  temperature  profile  with  the  following  shape:  T  - 
Tm  at  y  >  A/2;  T  -  T_m  at  y  <  -  A/2,  and  T  -  (Tw  +T_m)/2  in  between,  where  A 
is  the  minor  axis  of  the  elliptical  envelope  at  x  -  A/2.  This  shape  has 
been  measured  experimentally  by  Konrad  [32],  (see  also  Broadwell  and 
Breidenthal  [33],)  for  mixing  layer  flows  at  high  Reynolds  numbers.  This 
is,  to  our  knowledge,  the  first  time  it  has  been  computed  numerically. 

By  the  end  of  the  third  stage,  the  layer  cannot  absorb  any  more  energy 
and  a  relaxation  process  occurs,  during  which  the  kinetic  aind  thermal  energy 
are  fed  back  into  the  main  flow  streams.  This  reverse  action  is  accompanied 
by  the  fluid  leaving  the  core  and  moving  back  into  the  main  streams  at  a 
very  slow  rate. 

III. 3.  EFFECT  OF  MOLECULAR  DIFFUSION 

The  generation  of  large  temperature  gradients  within  the  core  as 
successive  instabilities  evolve  gives  rise  to  large  molecular  diffusion 
fluxes  which  act  to  smooth  out  some  of  these  gradients.  While  for  most 


where  6Q  is  the  core  radius  at  t  -  0  (for  more  discussion,  see  Leonard  [1], 
Ashurst  [34].)  The  cores  of  the  vortex  elements  and  of  the  temperature 
gradient  elements  become  different  with  time. 

Results  in  Fig.  23  show  the  temperature  profile  at  T  -  20  for  the  case 
of  X*  and  c  -  0.1 A*,  evaluated  for  a  -  0.0,  0.00001,  0.0001,  0.001,  0.01  and 

0.1.  Note  that  the  temperature  profiles  of  the  first  two  cases  are  almost 

identical,  indicating  that  the  effective  diffusivity  of  the  inviscid 
calculation  is  of  the  order  of  10-^.  in  the  last  case,  the  temperature 
profile  is  similar  to  the  case  of  pure  diffusion,  indicating  that  diffusion 
proceeds  at  a  rate  faster  than  the  instability.  It  is  also  noticed  that  fc 

moderate  values  of  a,  0.0001  <  a  <  0.01,  diffusion  only  affects  the  core  of 

the  eddies,  making  them  achieve  a  homogeneous  state  faster. 

Greengard  [35],  in  his  analysis  of  the  core-spreading  vortex  method  in 
which  a  fixed  number  of  elements  are  used  to  perform  the  convective 
transport  and  their  cores  are  expanded  to  account  for  the  effect  of 
diffusion,  showed  that  the  scheme  does  not  converge  to  the  correct  equation 
of  motion  except  when  the  flow  field  outside  the  region  |u>|  >  0  is  uniform. 
We  have  used  a  core  spreading  scheme  to  simulate  the  effect  of  diffusion  in 
the  energy  equation  with  two  modifications:  (1)  the  number  of  transport 
elements  which  discretize  the  gradient  field  is  increasing  with  time;  and 
(2)  a  is  kept  small.  Utilizing  an  increasing  number  of  elements  to  perform 
the  convective  transport  is  essential  since  it  is  important  to  determine  the 
gradient  field  accurately,  in  terms  of  the  location  and  strength  of  the 
elements,  before  the  aiffision  effect  can  be  added.  In  essence,  adding 
transport  elements  at  areas  of  high  strain  allows  the  computational  elements 
to  capture  all  the  vorticity,  and  temperature  gradient  carrying  fluid  at  all 
times,  even  after  the  vorticity  has  been  fragmented  by  the  action  of  the 
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strain  field.  Without  this  step,  strain  will  create  areas  which  are  void  of 
elements,  thus,  diffusion  cannot  be  represented. 

In  the  particular  application  of  a  shear  layer  at  high  Reynolds  number, 
the  flow  field  is  uniform  outside  the  area  where  |w|  >  0,  and  this  area 
expands  slowly  by  diffusion  if  a  <  1.  Limiting  the  simulations  to  small 
values  of  a:  (1)  reduces  the  errors  associated  with  the  fractional  step 
scheme  used  to  solve  Eq.  (24)  (Beale  and  Majda  [36]);  and  (2)  reduces  the 
errors  concomitant  with  convecting  an  element  with  the  velocity  evaluated  at 
its  center  while  its  core  radius  is  growing.  To  accommodate  this  growth, 
which  causes  the  spread  of  vorticity  in  the  direction  normal  to  the 
streamlines,  one  may  be  forced  to  add  elements  in  the  direction  normal  to 
the  maximum  principal  strain  direction,  and  then  redistribute  the  vorticity. 
It  is,  therefore,  clear  that  the  scheme  is  only  applicable  when  a  <  1  and 
for  short  time,  i.e.,  at  <  1.  If  these  two  conditions  are  not  satisfied, 
one  must  divide  each  element  whose  core  radius  is  larger  than  a  critical 
value  into  a  number  of  sepapate  elements  so  that  the  convective  transport 
can  be  performed  accurately.  Since  most  interest  in  shear  layer  flows  is  at 
high  Reynolds  number,  or  a  <  1,  and  within  the  short  time  of  development  of 
the  convective  instability,  we  feel  that  the  current  scheme  is  sufficient 
for  this  application. 

To  define  a  quantitative  measure  of  mixing  in  a  single  phase  fluid  with 
thermal  stratification,  we  observe  first  that  mixing  is  only  achieved  by 
molecular  diffusion.  Large  entrainment  fluxes  bring  the  unmixed  fluid 
layers  in  contact  along  a  larger  interface;  however,  molecular  diffusion 
across  this  interface  is  what  accomplishes  the  actual  mixing.  A  measure  of 
mixing  can  be  defined  as: 


M(t,a,e)  -  0;X_cD/“|T(x>t,a,c)  -  T(x,t,0, e)  |dx 

Note  that  M(t,0,e)  -  0,  while  M(t,a,0)  is  due  to  diffusion  only.  In  Fig.  24 
M( t, a, 0 .1 )  is  plotted  for  various  values  of  o  and  for  X*.  It  represents 
mixing  due  to  the  combined  action  of  entrainment  and  diffusion.  At  very 
small  values  of  a,  mixing  is  limited  by  the  amount  of  diffusion  across  the 
fluid  layers  which  have  been  entrained  into  the  eddy  core.  Since  for  these 
values  of  a  the  convective  transport  is  faster  than  the  diffusive  transport, 
mixing  increases  approximately  as  /a.  However,  as  a  increases,  and  at 
longer  times,  mixing  proceeds  at  a  slower  rate  since  it  becomes  bounded  by 
entrainment  of  unmixed  fluid  into  the  eddy  core  which  almost  ceases  by  the 
end  of  the  second  stage  of  rollup. 


IV  CONCLUSIONS 


In  this  work,  the  vortex  element  method  is  used  to  compute  both  the 
early  and  late  stages  of  development  of  an  inviscid  temporal  mixing  layer. 

In  this  method,  the  vorticity  is  initially  discretized  among  overlapping 
element  of  radially  symmetric  cores.  We  find  that  using  a  scheme  which 
depends  on  equating  the  vorticity  at  the  centers  of  the  elements  with  the 
accumulative  value  induced  by  all  elements  is  necessary  to  obtain  accurate 
results  for  initial  vorticity  discretization.  We  also  find  that  to  ensure 
the  accuracy  of  the  solution  for  short  times,  the  ratio  of  the 
core/separation  should  be  larger  that  one.  Very  large  cores  introduce  a 
strong  perturbation  in  the  vorticity  field,  while  smaller  cores  cause  a  fast 
deterioration  of  accuracy.  Using  fourth  order  Gaussian  cores  results  in 
better  accuracy  over  second  order  Gaussian  cores.  However,  we  feel  that  the 
improvement  in  accuracy  does  not  warrant  the  added  cost. 

As  time  proceeds,  the  distance  between  neighboring  elements  exceeds  its 
initial  values  due  to  the  generation  of  strong  stretch.  This  leads  to  the 
computation  of  inaccurate  velocities  and  is  manifested  by  the  irregular 
motion  of  the  vortex  elements.  To  overcome  this  problem,  the  vorticity  is 
constantly  redistributed  among  elements  inserted  along  the  principal 
direction  of  strain  to  capture  the  local  deformation  of  the  vorticity  field 
and  to  improve  the  resolution  of  the  calculations.  This  is  achieved  by  an 
insertion-and-interpolation  process,  which  is  applied  where  the  distance 
between  the  neighboring  centers  along  the  principal  direction  of  strain 
exceeds  a  threshold  value.  We  show,  using  solutions  for  a  shear  layer 
perturbed  at  different  wavelengths  and  amplitudes,  that  this  process  yields 
accurate  solutions  for  the  vorticity  distribution  at  long  times  and  after 
strong  strain  fields  have  caused  a  severe  distortion  of  the  streamlines. 


This  scheme  enables  one  to  accurately  compute  the  local  velocity  gradient 
which,  while  it  is  not  required  in  connection  with  vorticity  convection,  is 
necessary  for  the  accurate  evaluation  of  the  convection  of  a  passive  scalar. 

The  temperature  gradient,  distributed  over  transport  elements  which 
resemble  in  their  structures  the  vortex  elements,  are  used  to  compute  the 
temperature  distribution  as  the  rollup  evolves.  Contrary  to  vorticity, 
scalar  gradients  are  not  conserved  along  particale  paths,  thus,  the  strength 
of  these  transport  elements  is  changed  according  to  the  straining  and 
rotation  of  the  material  elements.  The  scheme  is  capable  of  capturing  very 
sharp  gradients  that  develop  within  the  core  since  the  elements  migrate 
towards  these  zones  by  convection.  The  multiplication  of  these  elements  via 
stretch,  which  inadvertently  mimics  the  physical  process  by  which  large 
scalar  gradients  are  gene reted,  provides  a  naturally  adaptive  grid  to 
compute  these  gradients.  By  expanding  the  cores  of  the  transport  elements, 
the  effect  of  small  diffusivities  can  be  simulated  as  a  small  perturbation 
to  the  convection  field.  Diffusion,  even  at  high  Peclet  number,  is 
responsible  for  generating  areas  of  uniform  temperature  inside  the  eddy 
since  it  acts  to  smooth  out  the  sharp  gradient  created  by  convection. 

The  application  of  vortex  methods  to  problems  in  which  the  no-slip 
boundary  condition  along  solid  walls  must  be  satisfied  can  be  accomplished 
using  the  random  vortex  method  (Chorin  [37],  and  Sethian  and  Ghoniem  [38].) 
In  this  method,  extra  vortex  elements  are  generated  along  the  solid  walls  to 
cancel  the  slip  velocity,  and  the  diffusion  of  vorticity  is  simulated  by  the 
random  walk  of  the  vortex  elements.  At  high  Reyolds  number,  a  strong  strain 
field  is  expected  to  cause  similar  problems  as  described  in  this  work,  i.e., 
areas  of  large  stretch  will  be  depleted  from  vortex  elements  and  accurate 
resolution  of  the  vorticity  field  may  be  lost  around  these  areas.  Extending 
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the  insertion-and-interpolation  scheme  described  in  this  work  to  the  random 
vortex  method  requires:  (1)  adding  a  third  fractional  step,  which  must  be 
performed  after  the  convection  and  before  the  diffusion  steps,  for  the 
redistribution  of  the  vorticity  field  among  elements  arranged  in  the 
direction  of  principal  strain;  and  (2)  computing  the  strain  field  at  the 
center  of  the  vortex  elements  in  a  Lagrangian  form  since,  due  to  random 
walk,  neighboring  vortex  elements  and  neighboring  material  elements  change 
as  time  evolves.  The  implementation  of  these  two  steps  must  be  preceeded  by 
careful  formulation,  and  will  require  lengthy  computation. 


KT. 

ft 
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FIGURE  CAPTION 


Figure  1.  The  growth  of  the  perturbation  amplitude  I  with  time  for  the  most 
unstable  case,  X*,  for  three  values  of  the  initial  perturbation 
e/X  -  0.001,  0.01  and  0.1,  showing  the  linear  range  and  the 
saturation  of  the  perturbation.  Each  curve  is  normalized  with 
respect  to  the  corresponding  value  of  I  at  t  -  0. 

Figure  2.  The  location  and  velocity  of  the  vortex  elements  during  the 

★ 

rollup  of  a  temporal  shear.  X  -  X  ,  with  e/X  -  0.001.  N(0)  - 

572,  h  -  0.3,  5  -  0.375,  and  At  -  0.1. 

Figure  3.  The  location  and  velocity  of  the  vortex  elements.  Wavelength  is 
X*,  and  e/X  -0.01.  At  t  -  0,  N  -  440,  h  -  0.33,  5-0.4  and  At 
-  0.1. 

* 

Figure  4.  The  location  and  velocity  of  the  vortex  elements  for  X  -  X  ,  and 
e/X  -  0.1.  All  the  numerical  parameters  are  the  same  as  in  Fig. 
3. 

Figure  5.  Schematic  diagram  showing  how  the  vorticity  is  redistributed 

among  three  elements  when  the  distance  between  two  neighboring 
elements  exceeds  a  pre-specified  value.  (xn,yn)  are  the 
coordinates  of  the  new  elements. 

Figure  6.  The  total  length  of  the  interface,  originally  at  y  -  0,  with  time 
for  the  cases  presented  in  Figs.  2,  3,  and  4,  normalized  with 
respect  to  its  length  at  t  -  0. 

Figure  7.  The  number  of  vortex  elements  used  to  represent  the  vorticity 

field  during  rollup  for  three  initial  perturbations,  normalized 
with  respect  to  the  corresponding  value  at  t  -  0. 

Figure  8.  (a)  The  total  kinetic  energy  of  the  perturbation  based  on  the 

2 

perturbation  velocity,  (U(x)-u(x,t) )  ,  and  (b)  The  total  kinetic 
energy  of  the  flow,  u^,  for  e/X  -  0.001,  0.01,  0.1. 


Figure  9.  The  evolution  of  the  vorticity  field  with  time,  compared  with  the 
experimental  results  of  Roberts  et  al.  [29]  for  the  spatial 
development  of  a  small  perturbation  of  a  shear  layer. 

Figure  10.  The  growth  of  the  perturbation  amplitude  for  X  -  10.5,  e/X  - 
0.01.  N( 0)  -  455,  h  -  0.3,  S  -  0.375  and  At  -  0.1. 

Figure  11.  The  location  and  velocity  of  the  vortex  elements  used  in  the 
calculations  of  the  case  shown  Fig.  10. 

Figure  12.  The  vorticity  field  for  X  -  2  X*.  e/X  -  0.01.  N(0)  -  540,  h  - 

0.44,  $  -  0.5,  At  -  0.1. 

* 

Figure  13.  The  vorticity  field  for  X  -  2  X  ,  e/X  -  0.1,  using  the  same 
numerical  parameters  as  in  Fig.  12. 

Figure  14.  The  vorticity  field  for  X  -  3  X*,  e/X  -  0.01.  N  (0)  -  818  and 

the  values  of  h,  6,  and  At  are  the  same  as  in  Fig.  12. 

★ 

Figure  15.  The  vorticity  field  for  X  -  3  X  ,  e/X  -  0.1,  using  the  same 
numerical  parameters  as  in  Fig  14. 

Figure  16.  The  location  and  velocity  of  the  vortex  elements  for  two 

★  ★  ★ 
perturbations,  X^  -  X  and  X^  *  2  X  ,  with  e  -  0.1  X  for  both 

perturbations.  N(0)  -  336,  h  -  0.55,  S  -  0.6  and  At  -  0.5.  A 

fourth  order  time  integration  scheme  is  used  to  transport  the 

elements. 

Figure  17.  The  total  amplitude  of  the  perturbation  of  the  case  in  Fig.  16. 

Figure  18.  The  temperature  distribution  across  the  layer  at  the  center  of 

* 

the  core,  for  X  and  e/X  -  0.001.  The  numerical  parameters  are 
the  same  as  in  Fig.  2. 


Figure  19.  The  temperature  distribution  across  the  layer  at  the  center  of 
the  domain  for  X  and  e/X  -  0.01.  The  numerical  parameters  are 
the  same  as  in  Fig.  3. 

Figure  20.  The  temperature  distribution  across  the  layer  at  the  center  of 
the  core  for  X  and  e/X  -  0.1.  The  numerical  parameters  are  the 
same  as  in  Fig.  4. 

Figure  21.  The  variation  of  the  logarithm  of  the  temperature  thickness  T 
with  time  for  the  cases  in  Figs.  2,  3,  and  4. 

Figure  22.  The  rollup  of  the  interface,  defined  by  the  layer  which  coinsides 
with  y  -  0  at  t  -  0  for  the  case  shown  in  Fig.  4. 

Figure  23.  The  effect  of  thermal  diffusion  on  the  temperature  distribution 
across  the  layer.  Temperature  is  shown  at  t  -  20  for  the  case 
shown  in  Fig.  4. 

Figure  24.  Total  mixing,  M(t,a,0.1)  due  to  the  combined  action  of 

entrainment  and  diffusion,  evaluated  for  different  values  of  a. 
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The  paper  on  "Numerical  simulation  of  a  reacting  shear  layer  using  the 
transport  element  method"  describes  the  formulation  of  the  vortex  element 
method  and  the  transport  element  method  for  a  chemically  reacting  flow  with 
finite  rate  of  heat  release  governed  by  Arrhenius  chemical  kinetics.  The 
application  of  the  method  to  a  ron-reacting  spatially  developing  shear  layer 
and  a  reacting  temorally  evolving  shear  layer  are  presented. 
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solution  of  tr.e  system  of  equations: 


u.  X.  i  *  th»-  v  :■  r  t . 


■ ;  es  expansion,  it  "<jr 


f  r. -°t  n-  c.;gr. 
,  a ' .  *  t r .  ■">  t  r  it.  ^ 


Fcr  ar.  inrompress ir.e  f  Et*,  Ec. 


is  c o r . s t a r. t  a . o r. g  par; 


*  X  x  X 


u  X  l  .  * 


where  x,  X  ,0;=*X  .  Tc  obtain  tn-.  v*.  Io:.t 

of  a  collect: on  of  vortex  elements  m  t r. e 
Eq.  v  1  w  ,  we  note  tnat  tr.e  stream  fur. otic 
single  vortex  element  is  Obtained  Zj  ir.it 
Eq.  (*♦  .  Using  po.ar  coordinates  t:  intt 


t.ms  equation 


vortex  element  p. 


dV.-'or  *  -<•.  r.-  d;  r,  w r.e-e  r  «  ..  r 

Tne  velocity  fie^d  of  a  single  e.emer.t 
radially  symmetric  a nr  u^=  -  o;.  cr. 

fielr  induced  by  a  distribution  of  fir 
vortex  elements,  cf  sr.apt  f.  ar.r  strer 

located  at  x  X  , t  IS: 

N 

u  ,x,t  .  '  :  • 


Ac?or2;ng  to  Eq.  ,‘E  ,  vortex  e : e-.--r.tr 
at  tr.e  lotaE  veiooity  at  tne;r  center?.  Ar 
progrerre-r,  tne  sirtnr.te  between  neigntt-'.n 
element?  increase*  sr.  tne  direction  of  max: 
?tra:n  rate  ?jcn  that  L i  >  r.(  wne-e  L>  i?  t 
distance  ir.  the  direction  of  maxsrcjir.  rfair. 
defines  a?  fix  ■  vfiu.fiy.  ]fiu'  anc  i  i.«  tr.e 
d i f fe-ence  operator  between  neighbo-:ng  ele 
Thi?  leas?  tc  a  deterioration  cf  the 
discretization  accuracy  since  accu-ate 
discretization  re^:-e?  that  {  >  i,.  Thu?, 
aigorjtnm  mu?t  be  uses  suer,  that  when  i>  s 
whe-e  f  •  1.E,  a  computational  element  i?  i 
at  the  midpoint  between  the  original  elemen 
circulation  of  the  ne„  element,  anc  tnat  o'. 
original  two  neighboring  elements,  1?  one  t 
the  sum  of  the  Circulation  cf  the  original 

e.eraents^  . 

For  com; resr i b.‘  ba-oc.mi'  f.ow,  E  ; . 
snow?  that  a  v  i  -St  •  E .  Moreover,  :•  «, 

where  A  i?  the  area,  wr.  ;.e  p  3A  *  ccr.sta.r 
Thu”,  the  o  i-Cu -at  i  or.  .?  constant  aior.g  a  ; 
path  --  Ke.vir.  theorem  --  anc  Ep?. 
usee  to  compute  the  evolution  of  the  v..rti. 
ve.ocity  f ie.s  provides  that  E3.  :}  is  use 
com;  at-  the  i-rct at  i  or.a  .  component  of  tr.»,  v 
due  tc  volumetric  expar.sior.,  a?  wi.,  be  Q-, 
tne  r.ex’  sect;;-..  When  ?;  1  V;  «  1,  tr.- 
c  lrc .,  1  at  icr.  cf  eacr.  vortex  e.emert  ».'!  b- 
eacn  tim*  step.  Esihg  the  definition  „f  tn 
ci-Ca.ation  ir.  Ep.  ft .  ,  wp  get 


ng  tc  tr.e  -  ow  Macn  nuna- 
Vp  c  =  -  VT^r,  wrue  71 


V. ,  :  ,  75 


S5  7T  nw 
5 how  ho« 

,  an  a:  i? 


ting  %r.~  velocity  of  tlu  vortex  eieme 
r.  craer  form . .  a .  Equat  lor.s  (1’}  ana 
tags airi  using  a  I  ourtr.  order  Hung-  ■ 

■.  method  w; tn  variable  fine  step  f  or 


Ar.oqr.s-  important  development  in  tn* 
at  i  or.  of  partio.B  netnoj?  to  reacting  f.cws 
■.re  fornvj.&lior.  of  tne  transport  element  method 
compute  tne  temperature  ant  species 
h-rntrat  I  or.  a  .  v.r.  r  or.s  ir.  a  Lagrangiah 

■  " '  .  1c  tr.:s  sent me,  trie  g-aaieno  of  tr.e 
el  -.r  fie.;  ;?  0  i  sere  t  i  ze  -d  into  a  number  of 
•it--  •-■-omer.ts  using  Fq.  witr.  u  rep. ace 2  by 

•  V?,  wr.ere  s  .?  a  generalize:  scalar,  oeir.g 
the.-  T  or  c.  Live  vortex  element?,  transport 
er-nt?  a r-.  a . s cr : cubed  w.nere  7?  >  2  aria  are 
■e:  wit:,  the  local  velocity  field  with  time, 
us ,  part.  .. e?  are  used  to  transport  sea.tr 
adients.  However,  contrary  to  vorticity,  scalar 
adients  are  not  conserved  along  particle  paths, 
d  should  be  modified  according  to  the  local 
r am  and  tilting  of  the  materia,  elements.  Tne 
tension  of  this  metnoa  to  reacting  flow  will 
|.i-t  changing  tne  gradient  transported  by  each 
emer.t  according  to  the  reaction  source  term  ir. 
s.  16,7,61  In  a  way  similar  to  changing  the 
-c..  ,t:ph  witr.  tne  no.n-barociir.ic  torque.  Thus, 
e  evolution  of  tne  enemies!  reaction  witr.  time 
..  oe  computed  ir,  a  Lagrangiar.  frame  of 
ferer.ee  as  the  interacting  species  flow,  in  tne 
-.owing,  we  describe  tr.e  conservative  form,  f 
e  transport  element  scheme  and  its  extension  to 
ive  Eqs.  it, 7, 5  . 

Initial. y,  tne  sea. or  gradient  g  .? 
scretized  according  to 


g  .  f ,  1  X  -  X 
-1  5  J  1 

1  =  1  J 


■-■re  .  x ,  c  and  h  have  tne  same  definitions  as 

fore,  an:  should  be  cnosen  to  satisfy  tne  sarn- 
quire-rents.  To  see  now  to  transport  the  scalar 
alien",  ir,  a  Lagrarigiar,  form,  we  start  by  trie 
n-ditfusice  case,  if  ?  is  a  passive,  non- 
iius.v  sca.ar,  tn*  conservat  ior.  eauat  lor.s  for  s 


=  car.  c-  se-.n  by  waiting  tn-  Giffe 


-  1  =51.  vu 


u  X 


g  fit.  -g  t  »  i  1  •.  t  *  1 1 ;  0  -  -  t  ,  Where  g  .  |g  . 
Trios,  tne  gradient  transported  by  an  element 
f  dial,..-'  a  parti  c.e  path  can  be  written  .»?:  g 

is,  i  t.  n  '  t )  r.  ~ ,  whiie  g  g  =  n  is  tr.e  ve  ■ 

norm  1.  to  il  '  ir.  a  chemical. y  reacting  flew, 
may  also  be  5s  :  .  For  a  graphical 

represer.tat  ion  of  tms  concept,  see  Fig.  1. 


i 

t  \ 

k  4  i 
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Figu-'e  Schematic  sweten  showing  the  evoiuti 
of  a  material  layer  separating  two  values  of  th 
conserved  scalar,  s,  ana  the  associated  sea. a- 
gradient,  g,  under  the  effect  of  stretch. 


Furthermore,  tne  fiux  initialized  by  Ed 
becomes: 


l  gjit)  -  Xj'.b 


rm  y  v  ’M* -j*  r 


x.  X:  ’  =Xr  i  l- 
a:;j  n  .  51  ; 


g  .  .  : 


n. j .—  _.  j  . 


to  u?Jd*.  g  t:  i:  L .: . 

ri:»tfve',l  tne  cxpr'e??  io:;  in  [  ; .  .. : 

t  T.  t?  Vjl’Jv  or'  5?  .  MO!”60  »'•?:' ,  i  n  « t  *ir -i  3 

.grpi*  i  r.g  e.  ^ .  . . _ _  i.  j  up n -1,  »  or.c  o  ar.  **<.iV 


:■.  mi,-.  :  *  ••  f-.’o.m  Jmce  7.=  =  g,  in.,?  ty 

mg  tr.v  grudir-.-.  7'?  *  V.g.  The  ?c.ut:on  oiy 

a?:  5  •  ;  V.g  j  ax,  wjwry  3  *  -*-u-  Xr.  r  :?  tr.r 

?:i.  m  f.ux,  V.g.  Integrating  by  par-.?,  pr.-  get? 
?  *  .'  g  .  V  j  :x.  '.'ring  Eq?.  '  -  -  >  ana  mi  far  g. 


g.  ’  ■  VX,  x-x.  'it. 


where  <  r)  «  /  r ’  fir'  :r',  a?  defined  before. 

If  pne  distance  between  .neighboring  element?  ir. 
tr.e  direction  of  prinoip.e  .‘train?  exceed?  a 
max infjir.  distance  Eh,  one  element  i?  inserted  half¬ 
way  betw-en  tne  two  element?  ana  tne  vi.ut  of  5. 
i?  adjusted  for  tne  three  element?.  A 
recombination  procedure  ear.  a.?o  be  implemented  to 
euro  tne  growt.n  in  the  number  of  computational 
e.e.T.e.nt?.  Tne  r.eoc  for  thi?  in?ertion  procedure 
i?  more  apparent  n ere  ?inee  the  magnitude  of  tne 
gradient  increase?  where  tne  strain  field  i?  hign; 
and  to  maintain  tne  acuracy,  more  element?  mu?t 
p-~;  'j 5 o "i  to  triTi5p jrp  i.c  gradient.. 

with  finite  diffusivity,  the  fir?t  term  on 
the  right  r,ar.d  ?iae  of  Eq?.  (6-8i  should  be 
?imj_ated  in  tne  ?olution.  in  gradient  form,  the 
conservation  equation  can  be  written  a?: 


~  -  -  g.Vu  -  i  v^g 


wn-.-re  i  i?  tne  mo-eojiar  diffusivity,  or  the 
i nverse  of  tne  Pec  let  njmoer.  At  hign  ?peed  flow, 

;  c; 

thi?  i?  typioaily  10"-10  .  To  ?oive  Eq.  (26 1 
u?ir;g  tne  ?cheme  tnat  we  have  developed  ?o  far, 

5i eier.er.t  g  mist  be  updated  according  to  tne 


diffusion  equation: 


li  =  a  V  g. 


progresses 
Peclet  nut 


-«s  of  tne  element?  will 


tot  a.  strength.  rtc 
Here  ?ir.Oe  We  l?e  i 


g  .  Vu  *  1  7'  g  -  1  crp  8 


where  <  i?  tne  number  of  ?pecie?.  1?  mg  tne 
definition?  of  g.  tne  strength  of  the  grad.  :.: 
mu?t  be  modified  according  to: 


d  ,  =  d w 

«  0?1  '  8*1  °?j  ?J1 


wnich  car.  be  written  a?  a  quadrature 
d ini :?>  at  *  dw'ds. 

Tne  a.gor ;t:.m  of  tne  tran?port  element  m*.  t:.. 
proceed?  a?  follow?:  l 1  .  update  the  location?  c: 
tne  eiement?  x.  according  to  tne  ve lOcity  at  tr.e i 

i 

center?  u?;r,g  Eq.  ! . 1  J :  (.  1  update  the  Vaiue?  of 
61  anc  n  either  according  to  the  integration  of 
Eq.  .....  or  oy  weeping  trace  of  tne  r.e  ignporing 
elements;  ,j,  update  tne  core  radii  of  different 
element?  according  to  the  corresponding  ? eo.et 
number  u?ing  Eq.  !2o);  and  (-  compute  the 
concentration?  of  ail  the  ?ca,ar?  u?ir.g  Eq.  (2-  • : 
and  ..p,  update  the  vnue  of  6?i  according  t.  Eq. 

(311.  in  mort  cases,  it  i?  po??ib.e  to  u?e  tne 
?ame  ?et  of  particle?  to  tranpport  element?  cf 
different  ?caiar?,  a?  well  a?  tne  vortex  e-em-mt; 
re?ulting  in  ?ub?tantiai  paving?  in  the  transpor 
?tep. 

Tne  tranpport  element?  generate  an  expan?:: 
field  a?  their  temperature  change?,  according  to 
Eq.  (j;.  The  velocity  field  a??ooiatea  with  tr.i 
expan? ion  „t  tne  .ow  Macn  number  limit  can  b- 
obtained  written  a?: 


Tne  total  velocity  produced  by  t 
field  i? : 

Vq.x.t.  -  I  V0«(X'Xi'i:O 

i»l  i  i 

wri'jre  VG  i.c  given  l>v  Ki.  ( G G  . 
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arranges  vertically 


discretize  tms 
patent is*  veioc 

aa 


vorticity  usir 
ity  component , 


^C,.  T.ne 
computes  Dy 


step-  *  n*.  -  *  a£)©rs-  -.■> 


'.£  two  source  flow?  o*.  x  ■  -*  ana  y  *  ♦;  ana  y 
*  -0  to  tne  velocity  fie.c  in  Eq.  ;2:  to  satisfy 
t.ne  boundary  cons  it  ion  at  x  -  0.  T.ne  nc-f.ow 
Dounda-y  condition  across  the  solid  wails  is 
implemented  Dy  using  conformal  mapping  arc  image 

,  *7 

vortices'.  At  the  downstream  sice  of  tne 
computational  window,  at  x  *  c,  vortex  elements 
are  aeieteo.  7n:s  induces  a  strong  perturoatior. 
wr.icr.  ensures  tnat  tne  rollup  and  first  pa;-ir.g 
will  always  taxe  place  within  tne  computational 
window.  Since  this  perturpation  is  not  applied  ir. 
an  organized  manner,  the  resulting  s.near  layer 
will  oe  considered  as  ar.  unforced  layer. 

Figure  2  shows  tne  location  and  velocity  of 
all  vortex  elements  used  in  tne  computations  for 
five  time  steps.  The  time  step  of  tne 
computations  is  it  *  C.i;.  Tne  picts  exhibit  a 
very  clear  and  accurate  portrait  of  the  rollup. 
During  rollup,  the  vorticity  within  the  shear 
layer  is  attracted  towards  the  center  of  a  large 
edoy,  entraining  fluid  from  both  sides,  and 
forming  what  appears  to  be  a  moving  focal  point  of 
a  spiral.  Between  neignbonng  large  eddies,  a 
zone  of  strong  strain  is  developing  where  the 
vorticity  is  depleted  and  the  gradients  are 
growing.  This  ’’braids"  zone  can  be  described  as  a 
moving  saddle  point  where  locally  the  fluid  flow 
experiences  a  separation  into  two  streams;  one 
moving  towards  tne  left  and  the  otner  moving 
towards  the  right  with  respect  to  the  saddle 
stagnation  point.  Downstream,  the  process  of 
rollup  continues  until  a  stronger  perturbation 
forces  two  neighboring  eddies  to  interact  ir.  a 
pairing  process.  it  is  important  to  stress  that 
the  algorithm  of  inserting  elements  as  the  strain 
field  develops  is  responsible  for  maintaining  the 
organization  of  the  calculation  for  a  long  time. 


Figure  2.  The  development  of  ia"ge  scale  vortex 
structures  in  a  spatia.  shear  layer.  Fact  •  t;r.t 
in  tne  figure  represents  a  vortex  element,  anc 
line  is  tne  velocity  vector. 


Quantitatively,  the  natural  frequency  cf 

shedding,  is  defined  as  f  *  U  /A,  where  2  * 

n  m  m 

(in*U2)/2,  and  A  is  the  wavelength  of  tne  largx 
eddy.  The  corresponding  Strouhai  number  is  S.= 

l/f  -  0.033.  This  is  the  same  value  as  the 
n 

frequency  of  the  most  unstable  mode  computed  free 
tne  linear  stability  theory  of  a  spatially 
developing  shear  layer  under  tne  cone  it  ions 
described  above.  Preliminary  results  for  tne 
growth  rate,  average  velocity  and 
turbulent  statistics  were  presented  in  study  of 
7 

Ghonierr.  and  Ng  for  the  forced  snear  layer. 

The  high  resolution  of  the  transport  element 
method  demands  tne  use  of  a  large  number  of 
transport  elements.  Moreover,  the  number  of 
elements  grows  rapidly  with  time  due  to  tne  severe 
stretch  produced  in  the  shear  layer.  This  maxes 
the  computation  of  a  wide  window  whicn  contains  a 
number  of  successive  eddies  expensive.  In  the 
next  section,  we  direct  attention  towards  a  mode, 
of  this  problem  that  requires  less  effort 
computationally  while  essentially  preserving  ai. 
the  physical  processes  involved  in  the  spatially 
developing  layer.  This  is  the  temporal  snear 
iayer  model  in  whicn  a  computational  window  that 
moves  at  the  average  speed  of  the  flow  is  impose: 
or,  a  single  wavelength  while  the  eddy  is  growing. 


v.  tempcrally-develc-im,  reaitin: 
she;:-  laye- 

Computational  results  snowing  tne  evolution 
of  a  large  easy  in  a  tempo-a .  sneer  layer  are 
presented  in  Figu-e  3.  In  tnis  case,  tne  boundar 
condition?  are  period::,  i.e..  w'x.y.t 
u'.x*A,y,t'.  and  u.x.y.t  -  u;x*»,y,t  ,  une-c  ;  ir. 
the  waveiengtr.  of  tne  perturbation.  Since 
detailed  analysis  of  trie  evolution  of  tnt 
temporal,  thermally  stratified  snear  layer  was 

presented  In  Gnor.iem  et  al.£  .  it  will  not  be 


is  Ciear. y  seer. 


presented  In  Gnor.iem  et  ai.  ,  it  will  not  be 
repeated  here.  Tne  qualitative  resemblance 
between  the  development  of  large  eddies  ir.  a 
spatial  and  a  temporal  snear  layer  is  clearly  seer, 
by  comparing  Figures  2  and  3-  Moreover,  tne 
snedamg  frequency,  i.e.  tne  frequency  of  the  most 
amplified  mode,  is  almost  tne  same  in  both  case. 
However,  the  growth  rate  of  tne  perturbation  is 
different  since  it  depends  on  tne  ve.ocity  ratio 
across  tne  layer;  a  parameter  tnat  does  not  appear 
in  tne  analysis  of  the  temporal  layer. 

In  tne  computation  of  tne  temporal  layer,  the 
window  is  limited  tc  or.e  wavelength  and  one  car. 
afford  to  use  more  elements  within  the  domain  to 
improve  tne  resolution.  One  car.  also  conduct. 
Inexpensive., ,  parametric  studies  on  tne  effect  of 
various  physical  parameters  that  appea-  in  tne 
mode.,  Eqs.kl-u,.  Tnus,  tne  temporal,  layer  wi.l 
be  used  as  a  model  for  the  spatial  layer  to  study 
turpulence-comoust  ion  interactions  ir.  shea-  flow. 


time  -  o.o: 


IME  «  11.01 


TIME  -  A.OC 


TIME  ■  16,00 


TIME  -  8.00 


TIME  -  20. 0C 


Figure  3-  Tne  development  of  a  large  eddy  in  a 
temporally  groming  shear  layer  at  the  same 
conditions  as  in  Figure  2. 


Since  tne  flow  in  unconfined,  tne  wavelength  i  is 
used  instead  of  H  to  non-dimens ional ize  the 
length . 

Tne  temperature  profile  across  the  midsection 
of  tne  eddy  is  exhibited  in  Figure  A.  The  roliup 
brings  fluid  from  one  side  to  the  opposite  side, 
while  stretch  increases  tne  gradient  across  each 


X  -  £.6  ,  TIME 


X  *  fc.t  .  TIM;  . 
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T 

X  -  6.6,  TIME  -  A. oo 


0 .  C  _  1.1 

X  -  6.6,  TIME  =  ’ 


X  *  6.6  ,  TIME  *  8.00 


X  -  6.6  ,  TIME 


Figure  Temperature  c istr iout i or,  a-rc.cc  f  - 
midsect  lor.  of  tne  large  eddy  shown  m  Figure  - 

layer.  Tne  proflies  show  that  after  tne 
relaxation  of  the  first  rollup,  a  secondary 
instability  develops  which  forces  tne  core  tnro.g 
another  turn,  tnus  creating  a  more  ragged 
temperature  distribution.  Moreover,  it  car  be 
seer,  that  tne  rollup  of  tne  shea-  layer  is  tr.e 
mecnanism  of  entrainment  tnat  leads  to  strong 
mixing  enhancement  as  the  two  fauics  diffuse 
across  tne  stretened  interface.  Since  ro.lu;  is 
associated  with  strong  stretch  tnat  reduces  tne 
tnicxness  of  tne  material  layers,  it  increases  ir- 
gradients  across  these  intertwining  layers,  thus 
ennancmg  the  diffusion  fluxes.  Quantitative.,,, 

the  rate  of  mixing  can  be  expressed  as  M  *  .  q  .  I 
da,  where  q  is  the  diffusion  flux,  n  is  tne  uni: 
vector  normal  to  the  material  surface,  and  da  is 
the  surface  area  element.  Since  for  twe 
dimensional  flow,  da  -  dl,  and  since  q  /  6.  - 

constant,  then  M  is  proportional  to  (£i)t.  Tne 

net  result  is  that  stretch  by  a  factor  ;  ennancer 
? 

mixing  by  a  factor  t  .  The  quadratic  rise  m 
mixing  during  rollup  will  have  a  significant 
effect  on  tne  rate  of  reaction. 

In  the  reacting  layer  calculations,  the  full 
system,  of  equations  is  integrated  using  partic.es 
that  transport  vortex  elements,  temperature 
gradient  elements,  reactant  and  product  gradients 
elements.  At  time  t  -  0,  the  vorticity  layer  and 
the  flame  front  coincide,  and  the  thicxnesse?  of 
the  vorticity  layer  as  well  as  the  temperature  an: 
species  concentration  within  the  layer  are  taKer 
to  be  equal.  A  small  sinusoidal  perturbation  wits 
amplitude  c  -  0.05  l  is  imposed  on  both 
distributions.  The  first  case  to  be  computed 
corresponds  to  the  following  set  of  parameters: 
pe"  Le"  1.  Aj.«  1,  Q  -  A  and  Ta-10,  and  n  -  l. 
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Figure  -  ?now«  tne  reacting  ?ne 
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However,  a?  rollup  ‘'tart?,  tne  folio 
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Figure  5.  Tne  development  of  a  .arge  eddy  ir.  a 
reacting  temporal  shea'-  layer  at  tne  sate 
conditions  as  ir.  Figure  • .  Tne  sc.io  .me  defines 
tne  flame  front. 


observed:  ( 1 )  a  swelling,  due  tc  tne  concomitant 
increase  in  the  rate  of  neat  release,  continues  as 
more  reactants  are  entrained  into  tne  burning 
core;  (2!  the  growth  of  tne  instaoility,  as 
measured  by  tne  ang.e  between  tne  major  axis  of 
the  elliptical  structure  and  the  main  stream 
direction,  is  encumbered  because  tne  volumetric 
expansion  causes  the  vorticity  intensity  to 
decrease  and  tne  eddy  to  become  weaxer  anc  less 
coherent;  and  (3)  tne  eddy  loses  its  symmetry  and 
becomes  eccentric  due  to  the  asymmetric  expansion 
and  the  generation  of  a  nor.-barocl  inic  torque 
associated  witn  heat  release.  As  more  of  tne 
initial  core  is  burnt,  the  fiuid  inside  the  eddy 
ceases  to  spin,  contrary  to  tne  nonreacting  case. 
Meanwniie  reactants  move  through  tne  side  to  enter 
the  reaction  region.  These  results  agree 
qualitatively  witr,  the  experimental  results  of 
Keller  ana  Daily5  on  tne  reacting  mixing  layer  at 
high  equivalence  ratios. 

On  tne  same  figure,  a  soiifl  line  is  plotted 
through  points  of  maximum  reaction  rate.  The  ’lne 
indicates  where  the  flame  front,  or  the  maximum 
heat  release  rate,  is  within  the  shear  layer. 

Below  this  line,  the  product  concentration 
approaches  unity  and  tne  temperature  reaches  T  . 

During  the  early  stages  of  rollup,  the  line  of 
maximum  reaction  rate  foiiows  one  of  the  material 
lines  closely,  i.e.,  the  growth  of  perturbation 


mrrt.y  Changes  the.  tope. eg,  cf  tr.-  f.a- 
At  later  stages,  this  l;n. ,  while  sea,-. 
ar.cir.e-  mate-io.  line,  form  a  tout  d/-, 
pro3iiv‘L5  wr.ere  tne  reactant r  a  -e  e:\tra* 
b.-h.  be  .ow  th  i «  lif.e,  me-*  p-Od-its 
co-e  almost  stops  its  rotation.  At  •*- 
of  burning  of  the  eddy,  tr.e  t»:  «.>.■«  ;• 
burr,  to  close  this  entry  way  and  tr.e  f. 

Oct  cf  tr.e  eddy  and  becomes  at,  t-rd.r.a',- 
flame. 

Tne  effect  of  heat  re. ease  or  tr.e 
of  the  easy,  wnior,  is  generated  by  the 
tne  sne a-  .aye*',  cat.  be  seer  from  the  te~' 
profi.es  across  tne  midsectior.  of  tr.e  wave 
Fig.re  t.  Since  tr. 
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Figure  6.  Temperature  distribution  across  the 
midsect  ion  of  the  reacting  large  eddy  shown  ir. 
Figure  5. 


core  of  tne  growing  eddy  from  the  right  side,  a  2- 
shaped  flame  is  formed.  At  the  initial  stages 
where  the  rate  of  entrainment  is  faster  than  the 
rate  of  burning,  the  flame  extends  deeper  into  the 
lower  stream,.  As  the  reactants  within  this  zone 
burn,  heat  is  released  within  the  core  of  the 
rotating  eddy,  causing  the  eddy  to  swell,  wniie 
maintaining  its  elliptical  shape.  The  non- 
barociinic  vorticity  generatea  around  tms  zone 
causes  the  observed  eccentricity  of  the  large 
edoy.  The  temperature  profiles  show  that  the 
higher  orde-  instabilities  observed  in  the 
nonreacting  case  are  suppressed  by  the  heat 
release,  and  that  the  core  of  the  eddy  stops  its 
rotation.  As  the  reactants  within  the  edoy  burn, 
tne  flame  leaves  the  structure  and  moves  into  the 
reactants.  This  results  in  the  formation  of  a 
temperature  profile  which  is  very  similar  to  the 
temperature  profile  at  t«0. 

To  study  the  effect  of  the  shear  iaye-  or.  the 
chemical  reaction,  we  compare  tne  rate  cf  burning 
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mass  of  products  generated  since  t  ■  C.  Figu-e  ■ 
depicts  the  tola,  mas?  Cjrr.,  »  ,  for  Doer  case , 


Figure  7.  Tots.  mas?  of  products  M  fornec  since 

P 

•  t  -  C,  in  the  reacting  shear  layer,  iabel.ee  a? 
RSL,  and  in  tne  laminar  flame,  laDelied  as  IF,  at 
tne  same  conditions,  and  tne  total  iengtn  of  tne 
flame  in  the  reacting  shear  layer  of  Figure  5. 

snowing  clearly  tne  different  stages  of  burning  in 
tne  reacting  shear  layer,  At  the  early  stages, 
during  the  linear  phase  of  development  where  the 
stretch  is  negligibly  smai.,  the  rate  of  burning 
is  linear  and  identical  to  that  of  a  lamina-' 
fiame.  As  the  layer  starts  to  roll  up,  the  area 
of  the  reaction  surface  increases  and  the  flame  is 
convoluted  around  tne  growing  eddy.  Tne  increase 
in  the  flame  area  due  to  convolution  is  nearly 
linear,  as  shown  in  Figure  7-  However,  as 
indicated  by  the  Figure  7  ,  the  products  form  at 
almost  a  quadratic  rate  du-'ing  this  stage.  Tms 
phenomenon  can  only  be  explained  by  recalling  that 
mixing  is  enhanced  as  the  square  of  the  stretch, 
and  that  in  this  case  of  fast  chemistry,  the  rate 
of  burning  is  governed  by  mixing. 

At  later  stages,  around  t  -  20,  the  amount  of 
products  forming  is  almost  nine  times  that  which 
forms  in  the  laminar  flame.  Due  to  flame 
convolution,  the  reacting  surface  area  has 
increased  by  three  folds.  Since  stretching  a 
layer  of  material  by  three  times  its  initial 
length  decreases  its  thickness  by  the  same  amount, 
and  the  fluxes  across  it  Increase  by  three  fold  as 
well.  This  augments  the  rate  of  mixing  by  nine 
times  over  the  non  stretching  case.  When  the 
chemical  reaction  is  fast,  as  in  this 

case,  the  material  mixed  reacts  and  the 
rate  of  reaction  Increases  by  the  same  rate  as  the 
mixing.  Tms  is  what  has  been  labelled  "mixing- 
controlled  reaction"  in  the  turbulent  combustion 
literature. 
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Figure  s.  Tne  temperature  7 

expansion  rate  e  a.ong  layer 
eddy. 


a t r q  i r.  ra:  e 
l r.  tr.c  rta 


To  further  analyze  the  res-.ts,  we  p.ot 
temperature  T,  the  strain  rale  s,  and  rat-.  :f 
expansion  e,  along  one  particular  layer  cf  f. 
within  the  reacting  eddy.  The  rate  cf  expar.= 
is  an  indication  of  the  rate  cf  product  forms 
i.e.,  tne  actual  reaction  rate.  Figures  8  ar. 


Figure  9.  The  temperature  T,  strain  rate  s,  anc 

expansion  rate  e  along  layer  8  in  the  reacting 
eddy . 
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Figure 
Figu-e  5. 


tr.e  reading  eddy  cf 


s.no*.  tnese  plot?  for  .ay  err  '  and  5,  reaped:  vtly , 
which  are  showr.  in  Figure  10.  Figures  6  arts  9 
snow  a  positive  corre-atior,  between  the 
temperature  and  the  strain  rate-,  eacr.  temperature 
pear  correspond?  to  a  maximum.  m  the  strain  rate 


VORTICITY  LOCATION 
STEP-  350  TIME-  17.57  LE-  1.0 


y  o- 


4 


SST-ik' _ --- 


Figure  11.  Tne  eady  in  a  reacting  temporal  shear 
layer  at  the  same  conditions  as  in  Figure  5  but 


t  -  17. 5'7. 


curve  and  vice  versa.  Moreover,  Figures  8  and 
9  exhibit  a  strong  correlation  between  the  local 
strain  rate  and  reaction  rate;  as  the  strain  rate 
decreases,  the  reaction  rate  decreases  ana  the 
temperature  drops.  These  results  indicate  that 
tne  rate  of  burning  and  temperature  are  positively 
correlated  with  the  strain  rate.  As  the  layer 
stretches,  the  diffusion  flux  of  tne  reactants 
into  the  flame  increases  and,  since  chemistry  is 
fast,  the  rate  of  burning  increases.  Under 
compression,  the  reactants  diffusion  fiux  is 
reduced  and  the  amount  of  burnt  mixture,  and  hence 
the  rate  of  expansion,  is  reduced. 


jt’CY'ed  =  i r.g  tr. 
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development  of  irv  _  e: 

snow:.  :r.  Figure  *  1 .  7 n*.  s 

reduced  ? : n c e  tnt  r^ie  cf  c 
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of  product? 


Figure  i; .  Tne  total  mass  of  pros u:ts  fc-r.v 


for  c  reacting  s.nea''  xayer 


C.5,  fil.,  an:  a  lar.ir.ar  fiane  at  tne  sam 


conditions, 


Figure  Tne  total,  mass  of  products  form? 


since  t  -  0  for  a  reacting  shear  layer 


0.28.  RSL,  and  a  laminar  fiame  at  the  same 


condition,  LF. 


parameters  shows  a  continuous  linear  rise  in  th 
mass  of  products.  If  tne  frequency  factor  is 


lowered  further  to  Af  -  0.25,  extinction  occur? 


earlier  at  around  t  -  10,  as  shown  in  Figure  '5 
Meanwhile,  the  corresponding  laminar  flame  shi« 
linear  rise  in  Mp .  Since  the  strain  rate 

increases  with  time,  the  extinction  phenomenon 
encountered  earlier  as  the  Damkohler  number  is 
reduced . 

To  explain  what  happens  around  extinction 


more  detail,  we  refer  to  plots  of  T,  s,  and  e, 
shown  for  layer  3  in  Figure  1 A  and  for  layer  s 
Figure  15,  both  for  Af  -  0.5.  Tne  plots  exhio; 

tne  variation  of  tne  three  parameters  at  t  ■ 
17.57.  The  geometry  of  the  two  layers  is  shown 


Figure  16.  Plots  for  e  show  that  the  expansion 
has  almost  stopped. 
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Figure  u,  Tne  teaperatu-e  T,  strain  rale  f,  ar.d 

expansion  rate  e  along  layer  5  ir.  the  reacting 
edcy  in  Figure  1 1 . 
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Figure  15.  The  temperature  T,  strain  rate  s,  ar 

expansion  rate  e  along  layer  5  in  the  reacting 
eddy  In  Figure  1  1  . 


Tnis  is  in  spite  of  the  fact  that  T  corresponds  to 
maximum  reaction  rate  along  most  of  the  layer. 
Contrary  to  the  case  of  Af  -  1.0,  the  values  of  T 

and  s  are  now  negatively  correlated,  l.e. 

temperature  maxima  correspond  to  minima  in  s  as 
shown  in  Figures  in  and  15  .  It  Is  also  observed 

that  s  and  e  are  negatively  correlated,  leading  to 
a  situation  in  which  strain  acts  to  extinguish  the 
flame.  This  Indicates  that  the  temperature  drops 


Figure  it. 
F  igu.-e  1 1  . 


Layers  3  end  5  m  the  reacting  ecc 


with  strain  due  to  the  increase  of  the  heat  f  1 
out  of  tne  layer  and  tne  reactants  f.ux  ir.t-  the 
layer .  Since  cnemical  reaction  is  s.ow,  tne  re---, 
loss  arc  the  gam  m  reactants  cannot  be 
compensated  by  chemical  heat  release,  .ead.ng  t: 
flame  blowout. 

vi.  conclusion; 

Advanced  numerical  methods  enable  one  t;:  i' 
integrate  elaborate  and  detailed  models,  wr.icr 
cannot  De  done  analytically,  sc  tnat  complex 
mechanisms  may  be  revealed  ano  analyzed;  arc  , _ 
provide  detailed  information  about  tne  flow  field 
which,  may  not  be  possible  using  traoiticna. 
experimental  techniques.  Computer  output,  rich  1- 
data,  offer?  a  challenge  in  how  to  extract 
valuable  information  about  the  phenomena  under 
investigation,  and  how  to  present  these 
information  in  compact  form.  Finding  out  tne 
appropriate  diagnostics  to  probe  computations, 
results  is  half  the  way  to  reaching  tne 
conclusions . 

In  this  article,  we  have  introduced  tne 
transport  element  method;  a  Lagrangian  particle 
scheme  based  on  the  discretization  of  the 
vorticity  and  the  gradients  of  the  scalars  intc 
finite  elements.  The  particles  move  along 
material  lines,  in  accordance  with  their  transpcr 
equations.  A?  strong  strains  develop  in  the 
dynamic  field,  the  finite  elements  may  change 
their  shape  or  configuration  to  accommodate  the 
distortion  which  is  produced  by  these  strain 
fields.  In  case  of  chemical  reaction:  (1!  the 
strength  of  the  elements,  i.e.  the  source 
strength,  changes  according  to  the  rate  of 
reaction;  and  (2)  the  chemical  heat  release 
induce?  volumetric  expansion  and  non-bar'ociinic 
vorticity  into  the  dynamic  field. 

The  simplest  model  which  can  Pe  proposed  to 
study  turbulence-combustion  interactions  contain 
five  parameters:  (1)  tne  Peclet  number  which 
defines  the  ratio  between  the  rate  of  convective 
and  diffusive  heating;  (2)  the  Lewis  number  wnicr. 
represent  tne  ratio  between  the  rate  of  heat  and 
mass  diffusion;  (3)  the  frequency  factor  whicn 
defines  the  ratio  between  the  rate  of  chemical 
reaction  and  mass  convection;  (A)  the  activation 
energy  of  the  reaction;  and  C 5 J  tne  enthalpy  of 
reaction.  The  outcome  of  these  interactions  can, 
thus,  be  presented  on  a  five  dimensions  spac^  on 


wr.icr.  one  can  xaer.tify  several  ? ut-d oma in?  Tor' 
b,,''r,:rg  enhancement,  f  laire  extincticr.,  f.a~ 

0=  :  : .  la:  ion? ,  etc.  To  accomplish  tr.  i?  goal, 
comp -tat  ion?  m-st  be  pe-forr.-o  f  o*-  a  r.a*.r;x  of 
para-eie'? .  The  compiled  data  Can  tner.  be  plotted 
on  tr.i«  spade .  Unce-  one  idealization  of  r.ign 
act  i  v  at  :or.  ene-g.  a  no  tr.  ir.  f.ame  structure  , 
results  cf  tr.e  asymptotic  ur.a.y?,?  car.  te  used  t . 
fi.I  ?o:-.  p  a  '  t  ?  cT  tr.:?  space  ar.o  snow  c.ne 

limiting  tren.c?  i~  '  i:  ’  . 

In  tr.i?  a-tic.e,  we  presented  result?  for  tr.e 

effect  of  Changing  tne  frequency  factor,  wrier. 

lead*  to  cnang.ng  tr.e  Dam>onj.er  nunoe",  at  fixed 

value?  of  tne  rent  of  tne  parameter?.  We  showed 

tr.ot  for  ?  ■  211,  L  ■  1,  T  *  1 1  arc  *  ■  - ,  at 
e  e  a 

A,*  '.C,  tne  ?tretcr.  a??ociatec  witr.  tne  ro-xup  of 

large  eddies  m  the  mixing  iayer  enhance?  tne 
combustion  proce??  by  a  factor  proportional  to  the 
equate  cf  tne  stretch.  7r.;«  l?  contrary  tc  tne 
a??umption  made  in  tne  conventional  w- i nx.ee 
lamina'  fiamr  theory  tr.at  the  inc'ease  m  tr.e  rate 
cf  Du-ning  if  linearly  proportional  t:  tne  streten 
factor.  Tr.  i?  quac'atic  increase  in  tr.e  teta.  rate 
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ABSTRACT 

The  vortex-scalar  element  method,  a  scheme 
which  utilises  vortex  elements  to  dise-etize  the 
region  of  high  vortitity  and  scalar  elements  to 
rep-esent  species  o-  temperatu-e  fields,  is 
utilized  in  the  numerical  simulations  of  a  two- 
d  :-er.s  i  oral  reacting  mixing  layer.  Computations 
a-e  performed  for  a  diffusion  flame  at  high 
Reynolds  and  Pe.let  numbers  without  resorting  to 
turbulence  models.  In  the  non-reacting  flow,  the 
mean  and  fluctuation  profiles  of  a  conserved 
scalar  show  good  ag-eemer.t  with  experimental 
measurements.  Results  for  the  reacting  flow 
indicate  that  for  temperature-independent 
kinetics,  the  chemical  reaction  begins  immediately 
downstream  of  the  splitter  plate  where  mixing 
starts.  Results  for  the  reacting  flow  with 
Arrhenius  kinetics  show  an  ignition  delay,  which 
depends  on  the  reactants  temperatu-e,  before 
significant  chemical  reaction  to  occurs.  Harmonic 
forcing  changes  the  structure  of  the  layer,  and 
concomitantly  the  rates  of  mixing  and  reaction.  In 
accordance  with  experimental  results.  Strong 
stretch  within  the  braids  in  the  non-equilibrium, 
kinetics  case  causes  local  flame  quenching  due  to 
the  temperature  drop  associated  with  the  large 
convective  fluxes. 


I  INTRODUCTION 

Turbulent  diffusion  flames  have  been  the 
subject  of  extensive  experimental  and  theoretical 
investigations  during  recent  years  (for  a  review, 
see  Bllger  [1]).  In  most  of  the  theoretical  work, 
turbulence  models  are  used  to  close  a  system  of 
averaged  transport  equations  which  describes  the 
statistical  behavior  of  the  aerothermodynamical 
variables.  Moment  methods  [2],  eddy  break-up  and 
mixing  controlled  models  [3],  flame  sheet 
approximation  [A],  assumed  probability  density 
function  (PDF)  shape  methods  [5],  solutions  based 
on  modelled  Joint  PDF  of  scalar  quantities  [6,7], 
and  based  on  modelled  Joint  PDF  of  scalar  and 
velocity  [8]  are  examples  In  which  turbulence 
modelling  have  been  used  for  the  closure  of 
equations  governing  the  statistical  quantities. 
Much  effort  has  gone  into  constructing  accurate 
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models  and  in  obtaining  results 
ag-eerents  wicr,  exper;ner  •. a  1  meas-.re---.ts. 
However,  in  complex  systems,  modelling  is 
difficult  because  of  our  lacx  of  knowledge  or.  tie- 
detailed  dynamics  of  the  flow.  Furthermore,  sine- 
most  of  the  interesting  dynarical  behavior  of  the 


flow  is  modelled  a 

pr 

iori,  such 

featjres  are 

not 

exhibited  from 

the 

results 

of  nur.e- 

i 

computations  based 

on 

turbulence 

models,  ar.d 

thus 

can  not  advance 

combustion. 

our 

understand 

:r,g  of  turbu 

lent 

The  progress 

in 

numerical 

methods  arc 

the 

availability  of  supercomputers  have  had  a  major 
Impact  on  tu-bulence  research.  Improved  accuracy 
of  the  numerics  and  increased  storage  arc 
computational  speed  have  made  it  possible  to  solve 
the  appropriate  transport  equations  governing 
turbulent  combustion  directly  without  the  need  re¬ 
modelling  for  some  limited  parameter  range.  Such 
model-free  "simulations,"  in  comparison  with 
calculations  utilizing  turbulence  models,  have  the 
advantage  that  the  physics  of  the  problem  is  not 
modelled  a  priori,  but  is  recovered  directly  from 
the  computed  results.  Their  results  can  be  used 
to  understand  many  important  mechanisms  of 
turbulent  transport  and  its  direct  influence  on 
chemical  reactions.  Furthermore,  since  the 
instantaneous  behavior  of  the  variables  are  known 
at  all  points  and  at  all  times,  accurate 
simulations  offer  a  good  method  of  probing  the 
flow  when  experimental  techniques  may  fail. 
Ultimately,  by  validating  the  results  of  the 
simulation  against  experimental  measurements,  ab 
Initio  predictions  will  be  a  reality. 

Numerical  methods  have  been  used  in  a  variety 
of  forms  for  the  simulation  of  turbulent  flows  in 
complex  configurations.  A  recent  survey  can  be 
found  In  review  articles  [9,10],  In  reacting 
flow,  three  approaches  are  used:  (1)  finite 
difference  methods,  (2)  spectral  methods;  and,  (3) 
vortex  methods.  In  the  first  approach,  the 
variables  are  defined  on  a  grid  and  the  transport 
equations  are  approximated  by  discretizing  the 
derivatives  on  the  grid  nodes.  Examples  of  this 
approach  can  be  found  in  the  work  of  Corcos  and 
Sherman  [11]  who  used  a  projection  method  to  study 
the  temporal  evolution  of  a  periodic  shear  layer, 
and  in  Grlnstein  et  al.  [12]  who  used  a  flux- 
corrected  transport  scheme  to  simulate  the 
development  of  coherent  structures  in  a  two- 
dimensional  spatially  evolving  shear  layer  and 
examined  their  effect  on  mixing. 

In  spectral  methods,  the  variables  are 
expanded  in  series  of  harmonic  functions  that 
satisfy  the  differential  equations  on  a  number  of 
collocation  points.  Riley  et  al.  [131  used  a 
pseudo-spectral  scheme  to  study  a  three 
dimensional  temporally-evolving  reacting  mixing 
layer  assuming  a  constant  reaction  rate,  constant 
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„■ : v :  et  a..  _  c.  use:  th-3  sat-  method  to  compute  a 
tw.  d o-.a 1  mixing  layer  with  ar.  A-rr,er.ius 
cr.e-.ical  reaction  and  constant  density  to  assess 
the  effects  of  large  coherent  structure  or.  the 
local  extinction  of  the  flare.  Extension  to 
spatially-growing  layers  was  initiated  fcy  Givi  and 
Jcu  [’tj  using  a  hytrid  pseudo-spectral  second 
orde-  finite  difference  scheme.  In  all  cases,  the 
Reynolds  number  was  kept  at  small  values.  Cl  1  Cel, 
limited  py  tne  g-id  resolution  and  the  r.umter  of 
ha-monic  modes. 

Ir.  the  th.i-d  approach,  vortex  methods  are 
used.  These  schemes  are  grid  free,  the  transport 
of  tr.e  variables  take  place  In  a  Lagranglan'  form, 
and  tne  solution  is  not  restricted  by  the  geometry 
of  tne  confinement.  Therefore  they  can  provide 
accurate  simulations  for  high  Reynolds  number, 
spatially  g-owing  flows.  Moreover,  vortex  methods 
optimize  the  computational  efforts  by  distributing 
computational  elements  a-ound  regions  of  high 
vorticity.  Tr.e  application  of  the  method  in  thin 
p-em.ixed  flame  calculations  with  a  finite  density 
jump  has  been  reported  by  Ghoniem  et  al.  [17]  and 
Sethian  [18],  among  others.  In  these 

calculations,  the  vortex  method  was  employed  to 
compute  the  flow  field,  and  the  dynamic  effect  of 
combustion  was  represented  by  the  propagation  of  a 
thin  Interface  at  the  laminar  burning  velocity 
acting  as  a  volumetric  source. 

Vortex  methods  were  also  used  in  simulating 
diffusion  flames  in  connection  with  a  finite- 
difference  approach  for  the  treatment  of  the 
scalar  variables.  Ashurst  and  Barr  [19]  used  the 
vortex  method  to  compute  the  hydrodynamic  field 
and  an  Eulerian  flux-corrected  transport  algorithm 
to  compute  the  diffusion  and  convection  of  a 
conserved  Shvab-Zeldovlch  scalar  approximating  the 
shape  and  convolution  of  the  flame  in  the  limit  of 
infinitely  fast  chemical  reaction.  Lin  and  Pratt 

[20]  used  the  random  vortex  method  to  simulate  the 
large-scale  motion  and  a  Monte-Carlo  method  to 
calculate  the  time-dependent  probability  density 
function  of  the  scalar  quantities  for  both  gaseous 
and  acqueous  mixing  layers.  The  PDF  transport 
equation,  however,  required  a  closure  model  for 
the  molecular  mixing  term. 

From  this  short  review,  it  is  clear  that 
numerical  simulations  have  played  an  important 
role  in  elucidating  the  physics  of  turbulent 
reacting  flows,  and  that  there  Is  a  continuing 
need  for  more  model-free  simulations  in  order  to 
explain  better  some  of  the  interesting  physical 
phenomena  that  have  been  observed  in  laboratory 
experiments. 

In  this  work,  we  extend  the  vortex  method  to 
study  non-premixed  chemical  reactions.  A  vortex- 
scalar  element  method  is  developed  to  treat  both 
the  hydrodynamic  and  the  scalar  field  in  a 
Lag-angian  sense.  The  fact  that  a  chemical 
reaction  Is  truly  a  Lagranglan  process,  i.e.,  It 
occurs  when  the  particles  (or  macroscopic 
elements)  Interact  as  they  flow,  motivate  the 
Implementation  of  Lagranglan  methods  fc- 
simulations  of  high  Reynolds  r.umbe-  r*j-i  :<v 
flow3.  The  method  Is  capable  --  ha-dli-v  i  - 
variety  of  Initial  and  boundary  •  edit  i  -  .>-  • 
not  limited  to  s,  :  >  f.  -  • 


and  Peel  et  nu-.ters.  In  this  paper,  we 
o-  tr.e  formulation  f  the  model  and  tr.e  numerical 
schemes,  and  p-esent  some  preliminary  validation 
studies  and  ir.te-pr  a  tat  ior.s  of  the  results. 

Ir.  Section  II,  the  geometrical  corf  i  g  „-at  i  c  r. 
of  a  spatially  evolving  mixing  lave-  is  presented, 
and  the  formulation  of  the  prct.er  and  cf  tr.e 
scheme  are  described.  Results  of  some-  sample 
calculations  are  given  in  Section  Hi. 
Computations  of  a  non-reacting  mixing  layer  is 
performed  first  in  orde-  to  check  on  the  accuracy 
of  the  method  by  comparing  Its  results  with 
experimental  measurements  at  the  same  conditions. 
Preliminary  results  of  a  simulation  of  a  reacting 
mixing  layer  in  which  the  two  reactants  are 
introduced  in  different  streams  are  presented 


next.  Both  constant  rate  kinetics  and  tempe-ature 
dependent  kinetics  are  considered.  Ir.  both  cases, 
the  influence  of  the  coherent  structures  on  tr.e 
finite  rate  chemistry  is  assessed  and  in  the 
second  case,  the  non-equilibrium  effects  ir.  the 
reaction  rate  are  examined.  In  the  constant  rate 
kinetics  calculations,  the  influence  of  harmonic 
forcing  at  the  inlet  of  the  mixing  laye-  is 
investigated.  This  study  was  motivated  by  recent 
experimental  observations  of  Roberts  and  RosV.: 

[21]  and  numerical  computations  of  Ghoniem.  and  Ng 

[22] .  The  paper  is  concluded  ir.  Section  IV  wit!",  a 
summary  of  our  new  results  and  suggestions  for 
future  developments. 


II  FORMULATION  AND  NUMERICAL  SCHEME 

A  two-dimensional,  confined,  plana 
layer  is  considered.  A  schematic  diag-a- 
flow  field  is  shown  in  Fig.  1.  Two  i 

unmixed  reactants,  fuel  F  and  oxidant 
present  at  small  concentrations  in  the 
speed  stream  and  bottom  low  speed 

respectively.  We  make  the  following  assu 
(1 )  tne  neat  release  is  low  so  that  its  e 
the  dynamics  of  the  flow  is  negligible; 
Mach  number  is  small;  (3)  the  free 
concentrations  of  F  and  0  are  equal  and  c 
(k)  the  molecular  diffusivities  are  e 
constant;  (5)  the  viscosity  is  the  s--- 
streams;  and  (6)  the  chemical  react ; 
and  0  is  single  step,  Irreversible. 
order.  The  density  is,  therefore,  • - • 
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the  transport  equations  of  the  hydrodynamic  field 
and  the  scalar  —  temperature  or  species  —  fields 
are  decoupled.  The  equations  governing  this 
system  are: 


F  +  0  — -  P 


is  unity,  another  conserved  scalar  can  be 
introduced,  BpT-  cp-  T/Q,  and  the  solution  of  Eqs. 

(6)  and  (7)  for  c  and  B  will  determine  the 
P 

behavior  of  all  the  scalar  quantities,  cp,  cQ, 
cp,  and  T. 
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II.  1  THE  VORTEX  SCHEME 

In  the  vortex  method,  the  vorticity  field  is 
represented  by  a  finite  number  of  vortex  elements 
of  finite  cores: 

u( x, t)  -  I  jy 62  f ( x  -  Xi)  (8) 

where  r,  -  J  w  dA,  is  the  circulation  of  a  vortex 


element  and  6  is  the  core 
center  of  the  element,  f 
distribution  associated  wi 
the  core  function  (Chorin 
Beale  and  Majda  [25].) 
obtained  by  solving  Eq. 
vorticity  distribution: 


radiu3,  while  Xj  is  the 
represents  the  vorticity 
th  a  vortex  element,  or 
[23]  and  Hald  [2H] ,  and 
The  velocity  field  is 
(2)  U3ing  the  discrete 


where  P  indicates  products  and  W  *  cpcQ  exp(- 

Ta/T)  is  the  reaction  rate,  written  in  terms  of 
the  rate  of  generation  of  products  per  unit  mass, 
u  -  (u,v)  is  the  velocity,  x  -  (x,y)  and  x,  y  are 
the  streamwlse  and  cross  stream  directions, 
respectively,  t  is  time,  *  is  the  stream  function 
defined  such  that  u  -  3i|i/3y  and  v  -  -  3*/3x,  u  -  ? 
x  u  is  the  vorticity,  e  is  the  concentration  per 
unit  mass,  T  is  temperature.  V-( 3/3x, 3/3y) ,  and 

V2-32/3x2*32/3y2.  Variables  are  non 

dimensionalized  with  respect  to  the  appropriate 
combination  of  the  total  shear  aU  -  1)1-  'J2,  the 
channel  height  H,  the  free  stream  concentration  of 
F,  cpo,  the  free  stream  temperature  at  x  •  0,  To. 

In  Eq.  (5),  J  -  F  or  0  for  fuel  and  oxidizer, 
respectively.  Re  -  nU  H/v  is  the  Reynolds  number, 
where  v  is  the  kinematic  viscosity.  The  reaction 
rate  constant  k  -  A  exp(-Ta/T)  where  A  is  the 
frequency  factor,  and  Ta  is  the  activation  energy, 
non-dimen3ionalized  with  respect  to  (RTo),  R  being 
the  gas  constant.  Q  is  the  enthalpy  of  reaction, 
non-dlmensionalized  with  respect  to  CpTo,  where  Cp 
is  the  specific  heat  at  constant  pressure.  Pe  * 
aU  H/a  is  the  Peclet  number,  where  a  Is  the 
thermal  diffuslvity.  Da  -  A  cpo  H/aU  is  the  first 

Damkohler  number.  Le-  a  /D  is  the  Lewis  number. 

Since  Eqs.  (if).  (5)  and  (6)  are  similar, 

there  is  no  need  to  solve  them  all  if  the  scalar 
concentrations  cp,  c^  and  cp  are  normalized  in 

such  a  way  that  their  initial  and  boundary 
conditions  are  identical.  This  is  accomplished  by 
the  use  of  Shvab-Zeldovlch  transformation  [1], 
Introducing  conserved  scalars  8PD  •  c  ♦  c_,  and 

p  P  r  P 

60P*  (co"  V*  He  eet! 


♦  u  .  V8. 


-J_  728 

Pe  Le  8J 


for  J  •  FP  or  OP.  Since  Bpp  and  8 
initial  and  boundary  conditions,  B, 


Since  8pp  and  6Qp  have  the  same 


- ,  - — ,  „pp  „op 

The  finite  rate  kinetics  effects  can  be  taken  into 
account  by  considering  the  transport  equation  for 
the  product  of  chemical  reaction,  Eq.<6),  and  Eq. 
(7)  for  a  conserved  scalar.  If  the  Lewis  number 


u  ■  I  Tj  Ktx-jj)  k(x-Xj)  ♦  up  (9) 

2 

where  IC(x)  ■  -(y,-x)/r  13  the  kernel  of  the 

Poisson  equation,  <(x)  -  /  r  f(r)  dr  is  the 
circulation  within  r,  and  r  -  |x|.  up  is  an 

irrotational  velocity  field  added  to  satisfy  the 

2 

potential  boundary  condition;  up  »  V$  where  V  $  - 

0  and  u.n  -  0  on  solid  boundaries  while  u.n  -  U  at 
the  inlet,  n  is  the  normal  unit  vector.  For  the 
confined  shear  layer,  the  boundary  condition  at  x 
-  0  is:  u  -  UI  for  y  >  0  and  u  «  U2  at  y  <  0, 
while  y  -  0  is  a  vortex  sheet  of  strength  flU-  U1- 
U2. 

In  this  work,  we  use  Ranklne  vortex  elements, 
i.e.,  the  vorticity  of  an  element  is  constant 
within  the  core  and  zero  outside,  f(r)  -  1/x  for  r 
S  6  and  f(r)  -  0  for  r  >  6.  Correspondingly,  x(r) 
2 

•  r  /2x  for  r  S  6  and  k  •  1  for  r  >  6.  Moreover, 
the  potential  velocity  field  is  obtained  by 
conformal  transformation.  Thus,  the  physical 
plane  is  mapped  onto  the  upper  half  plane  and 
image  vortices  are  used  to  satisfy  the  potential 
boundary  conditions.  The  form  of  the  mapping 
function  for  the  confined  shear  layer  is  given  by 
Ghoniem  and  Ng  [22]. 

The  motion  of  the  vortex  elements  must  be 
constructed  such  that  the  vorticity  field 
satisfies  Eq.  (3).  This  is  accomplished  by 
solving  this  equation  in  two  fractional  steps: 

Convection:  ♦  u.Vw  -  0  (10) 

Diffusion:  |f  *  ^  ?2“  (11) 

In  the  first  step,  the  convective  transport 
of  vorticity  is  Implemented  in  terms  of  the 
Lagrangian  displacement  of  the  vortex  elements 
using  the  current  velocity  field  computed  from  Eq. 
(9).  In  the  second  step,  the  solution  of  the 
diffusion  equation  is  simulated  stochastically  by 
the  random  walk  displacement  of  the  vortex 
elements  according  to  the  appropriate  population. 
Thus: 

Xj(t*at)  -  xt(t)  ♦  I  u(xlk)at  *  nj  (12) 


Reaction 


(17) 


for  i  •  1,2,...,N,  where  I  is  a  k-th  order  time- 
integration  scheme  and  xij  is  a  two  dimensional 
Gaussian 

random  variable  with  zero  mean  and  standard 
variation  /2at/Re.  For  more  details,  see  Ghoniem 
and  Ng  [22],  Ghoniem  and  Gagnon  [26]. 

The  no-slip  boundary  condition  at  the  walls 
is  satisfied  by  generating  new  vortex  elements  to 
cancel  the  induced  velocity  by  the  vorticity 
field.  Here,  we  generate  vorticity  only  at  the 
point  of  separation,  i.e.  at  the  tip  of  the 
splitter  plate  since  the  growth  of  the  boundary 
layers  along  the  channel  walls  at  these  high 
Reynolds  numbers  is  small.  At  each  time  step,  the 
new  vorticity  sf  -  -aU  Um  at,  where  Um  -  (U1  + 
U2)/2,  is  consigned  to  No  elements  of  strength 
ar/No  and  added  to  the  field  at  points  ax  -  Um/No 
apart  downstream  of  x  ■  0. 

The  effect  of  the  numerical  parameters  on  the 
accuracy  of  the  results  was  investigated  by 
Ghoniem  and  Ng  [22].  Their  results  emphasized  the 
importance  of  using  a  high  order  time-integration 
scheme  with  k-2  to  avoid  excessive  numerical 
diffusion  in  the  vorticity  field.  The  value  of  No 
-  6  wa3  also  found  to  be  appropriate  in  order  to 
obtain  well-defined  eddy  structures  after  the 
rollup  and  the  first  two  pairings.  The  second 
pairing  is  accomplished  within  the  domain  of  0  £  x 
£  6,  therefore  the  computational  domain  was 

limited  to  Xmax  -  6.  Downstream  of  Xmax,  the 

vorticity  was  deleted.  Varying  Xmax  showed  that 
the  effect  of  deleting  the  vortex  elements 
propagates  about  one  channel  height  upstream, 
hence  the  results  are  accurate  only  for  0  £  x  £5. 

II. 2  THE  SCALAR  ELEMENT  METHOD 

In  this  scheme,  which  is  a  two  dimensional 
extension  of  the  random  element  method  of  Ghoniem 
and  Oppenheim  [27],  the  scalar  field  is 
represented  by  a  set  of  elements  each  carrying  a 
finite  amount  of  the  scalar  field: 

s(x,t)  ■  I  Sj  6(x~xl)  03) 

where  s  is  a  scalar  field,  being  the  temperature 
or  species 

concentration,  s^  is  the  strength  of  an  element, 

defined  as  the  amount  of  scalar  carried  by  this 
element  and  s(.)  is  the  Dirac  delta  function.  Sj 

•  1/SA  /  s(x,t)  dA,  where  6A  -  6x6y,  and  Sx  and  Sy 
are  the  distances  between  the  centers  of 
neighboring  elements  in  the  streamwise  and  cross 
stream  directions,  respectively,  and  is  the 

center  of  the  element.  If  s  is  an  active  scalar, 
its  transport  is  governed  by: 

If  ♦  u  •  ?s  ■  4-  V2s  ♦  W  (10) 


Convective  and  diffusive  transport  are  taken 
into  account  in  a  similar  way  as  in  the  vortex 
method,  i.e.  by  the  Lagrangian  motion  of  the 
scalar  elements  using  the  velocity  field  u, 
computed  using  Eq.  (9),  and  the  random  walk 
displacement  of  the  elements  using  a  set  of 
Gaussian  random  variables  with  zero  mean  and 
standard  deviation  /2at/Se  (Ghoniem  and  Sherman 
[28].)  If  j.  is  the  center  of  the  element  1, 
then, 

X^t+at)  -  x^t)  +  E  u(x1(<)  *  ''i  (15) 

Chemical  reaction  changes  the  amount  of 
reactants  carried  by  the  element  according  to  the 
integration  of  Eq.  (17), 

s^t+at)  -  s^t)  +  W  at  (19) 

However,  the  reaction  occurs  only  when  the 

elements  are  close  enough  for  molecular  mixing  to 
affect  their  composition.  Therefore,  at  every 
time  step,  the  distance  between  the  centers  of 
each  two  elements  of  F  and  0,  axj.-l  Xj^Xj  I,  is 

computed.  If  ax^S  6D>  where  6D-0(1//Se)  is 

the  diffusion  length  scale,  the  composition  of 
each  of  the  two  elements  changes  according  to  Eq. 
(19).  The  initial  distance  between  neighboring 
elements  must  be  small  enough  to  allow  enough- 
interactions  between  the  elements.  Thi3  limits 
the  maximum  value  of  the  Peclet  number  that  can  be 
economically  used  in  the  computations  to  0(1000). 

The  scheme,  while  providing  an  approximate 
solution  of  Eq.  (12)  in  a  stochastic  sense,  mimics 
closely  the  actual  physics  of  the  reaction 
process.  This  is  achieved  by  using  the  lagrangian 
formulation  of  the  transport  equations  and  dealing 
with  the  chemical  production  terms  in  Individual 
particles. 


Ill  RESULTS  AND  DISCUSSION 

The  computer  code,  developed  by  Ghoniem  and 
Ng  [22]  for  vortex  simulation  of  a  non-reacting 
shear  layer,  was  vectorized  in  order  to  take 
advantage  of  the  computational  capability  of  a 
CRAY-XMP.  The  scheme,  being  explicit  in  time  and 
requiring  mostly  non  recursive  computations,  can 
utilize  this  capability  efficiently.  The  dynamics 
of  the  non  reacting  layer  was  investigated  in  ’ 
detail  in  the  work  of  Ghoniem  and  Ng  [22].  Here 
we  concentrate  on  results  pertaining  to  mixing  and 
to  a  chemica.ly-reacting  layer. 


where  Se  is  the  ratio  between  the  diffusive  and 
convective  time  scales  of  transport  of  s,  Se  -  Pe 
for  s  -  T,  and  Se-  Pe  Le  if  s-c.  In  the  scalar 
element  method,  this  equation  is  solved  in  three 
fractional  steps: 

3s 

Convection  -rr  *  u  •  Vs  •  0  (15) 


III .  1  NON-REACTING  MIXING  LAYER 

Results  of  a  typical  simulation,  presented  in 
terms  of  the  velocity  and  location  of  all  vortex 
elements  used  in  the  computations,  are  shown  in 
Figs.  2,  3  and  9  for  the  cases  of  Re  -  29000,  Re  - 
9000,  and  Re  -  1000,  respectively.  Each  vortex 
element  is  depicted  by  a  point,  while  its  velocity 
relative  to  the  mean  velocity  is  represented  by  a 
line  vector  starting  at  the  center  of  the  vortex 
element.  The  velocity  ratio  across  the  layer  at 
the  inlet  is  U2/U1  -  1/3. 


Diffusion 


(16) 


TINE  •  it.  to¬ 


ll*  •  ft.  00. 


«b  •  m  •  mom 

TINE  •  M»«0. 


«  •  ».  mi.  m  •  mom 

V1K  *00.00. 


0*o »n.  m  -  looM 


in  -  a  mo.  is  •  tooM 

TIM  •  n.ao 


MB  •  lltn.  IB  •  100M 
TIK  •  Ml  00 

Figure  3-  Vorticity  field  at  Re  *  10000, 
U2/U1  -  1/3. 


Results  show  the  formation  of  large  vortex 
eddies  by  the  rollup  of  the  vorticity  layer  that 
emanates  at  the  splitter  plate,  and  the  subsequent 
pairings  of  these  eddies  Into  larger  structures. 
The  rollup  of  the  shear  layer  was  Investigated  In 
Ghonlem  and  Ng  [22]  by  analyzing  results  at  a  wide 
range  of  the  Reynolds  number  and  at  different 
boundary  conditions.  Their  analysis  show  that: 
(1)  the  rollup  Is  due  to  the  growth  of 
perturbations  by  the  Kelvln-Kelmholtz  Instability 
mechanism,  and  the  shedding  frequency  corresponds 
to  the  most  unstable  frequency  predicted  from  the 
linear  stability  analysis  of  a  spatially  growing 
layer;  (2)  pairing,  which  Is  associated  with  the 
local  subharmonic  perturbations,  results  In  a 
step-wise  Increase  In  the  size  of  the  vorticity 
layer  as  two  eddies  merge;  (3)  The  two  sources  of 
the  subharmonle  perturbations  are  the  downward 
motion  of  the  layer  and  the  monotonic  growth  in 
the  size  of  the  eddies  downstream;  ( H )  the 
Intrinsic  dynamics  of  the  instability  is  not 
strongly  affected  by  the  value  of  the  Reynolds 
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Figure  4.  Vorticity  field  at  Re  *  4000, 
U2/U1  -  1/3- 


number,  except  that  at  the  low  Reynolds  number  the 
eddies  are  slightly  larger  due  to  the  dispersion 
of  vorticity  by  diffusion;  and  (5)  the  computed 
velocity  statistics  show  good  agreements  with 
experimental  data,  indicating  that  the  fundamental 
mechanisms  of  the  shear  layer  are  two  dimensional 
and,  hence,  the  numerical  scheme  is  capable  of 
predicting  the  large  scale  features  accurately. 

To  study  entrainment,  a  passive  conserved 
scalar  with  a  normalized  concentration  value  equal 
to  zero  In  the  high  speed  stream  and  equal  to  one 
in  the  low  speed  side  Is  Introduced  at  the  inlet 
section.  At  each  time  step,  19  elements  are 
introduced  in  each  stream.  The  initial  distance 
between  two  neighboring  elements  in  the  cross 
stream  direction  Is  taken  as  6y  •  0.021.  The  time 
step  at  •  0.1,  thus  the  distance  between  the 

elements  in  the  streamwise  direction  Is  6x  -  0.05 
on  the  average.  Since  diffusion  Is  more  critical 
In  the  cross  stream  direction,  6y  Is  chosen  to  be 
smaller  than  6x.  A  case  with  6y  »  0.016,  using  25 
elements  In  each  stream  was  computed,  showing  no 
significant  change  In  the  overall  behavior. 

Figures  5,  6  and  7  are  obtained  for  Reynolds 
number,  Peclet  number  and  velocity  ratio  10,000, 
4,000  and  1/2,  respectively.  Figures  5  and  6  show 
the  velocity  and  location  of  all  the  vortex  and 
scalar  elements  respectively,  while  Fig.  7 
exhibits  the  strength  of  each  of  the  scalar 
elements  at  the  non-dimensional  times  of  t  •  28, 
29  and  30.  In  Fig.  6,  the  dots  represent  the 
fluid  from  the  high  speed  side  with  normalized 
concentration  c  •  0,  and  the  open  circles 
represent  the  fluid  from  the  low  speed  side  with  c 
•  1.  This  figure  indicates  that  the  rollup  of  the 
vortices  and  their  subsequent  pairing  entrains 
fluid  from  both  sides  of  the  free  streams  Into  the 
cores  of  the  vorticity  layer,  which  results  in  the 
enhancement  of  mixing  between  the  two  streams. 
Entrainment  asymmetry  Is  observed  as  more  fluid 
from  the  high  speed  side  Is  present  In  the  low 
speed  side  than  the  opposite  (Koochesfhani  [29]). 

The  Instantaneous  profiles  of  the 

concentration  field  are  averaged  over  a  long-time 

period  and  the  statistical  values  are  compared 

with  experimental  data  in  Figs.  8  and  9.  Figure  8 

shows  the  mean  value  of  the  concentration,  c  ,  as 

in 
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a  function  of  (y-y_)/(x-x  ),  where  y  is  measured 
o  o  o 

at  c  -  0.5  and  x_  is  the  virtual 
m  o 

origin  of  the  mixing  layer  based  on  the  mean 
concentration  profile  (in  the  calculation,  x  - 
0).  In  this  figure,  the  solid  line  is  the 
computed  mean  concentration  at  x  «  4  and  the  data 
points  are  obtained  from  recent  experimental 
measurements  by  Masutani  and  Bowman  [30]  for  a 
dilute  non-reacting  mixing  layer  with  the  same 
velocity  ratio.  Figure  9  shows  a  comparison 
between  the  computed  and  measured  mean 

—  2  —  2 

fluctuations  of  the  concentration,  c’  -  (c-c  )  . 

Ul 

It  is  evident  from  the  two 

figures  that  both  the  mean  and  the  second  moment 
of  the  conserved  scalar  across  the  width  of  the 
shear  layer  are  accurately  predicted  by  our 
computations. 
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Figure  6.  Scalar  a  velocity  field  at  Re  -  10000, 
U2/U1  -  0.5. 
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Figure  9*  Normalised  rms  concentration  profile  as 
a  function  of  the  cross-stream 
coordinate . 


•  a  mi.  he  •  imm 


He  note  cnat  the  results  in  Figs.  8  and  9  are 
in  better  agreement  with  experimental  data  than 
those  previously  predicted  by  Givi  et  al.  [31]. 
In  these  calculations,  a  k-e  turbulence  model  and 
a  gradient  diffusion  model  for  turbulent  transport 
of  the  scalar  mean,  moment  and  probability  density 
function  was  utilized.  In  the  k-e  calculations, 
the  concentration  fluctuations  exhibit  a  fairly 
smooth  bell-shaped  profile  with  a  much  less  clear 
double  "hump"  in  the  middle  region,  Indicating 
poor  agreement  near  the  high  speed  stream.  The 
present  calculations  show  the  two  local  maxima  in 
the  fluctuation  profiles  that  correspond  to  the 
location  where  the  gradient  of  the  mean  value  is 
highest.  The  same  behavior  is  observed  by  the 
experimental  results  of  Masutanl  and  Bowman  [30] 
and  Batt  [32].  It  is  clear  that,  in  accordance 
with  the  findings  of  Broadwell  and  Brledenthal 
[33],  the  intermittency  caused  by  the  large 
coherent  structures  contributes  greatly  to  the 
statistics  of  turbulent  flows. 

III. 2  REACTING  MIXING  LAYER 

In  the  calculation  of  a  reacting  mixing 
layer,  two  reactants  F  and  0  are  Introduced  on 
both  sides  of  the  splitter  plate.  At  x  -  0,  for  y 
>  0,  Cp“  1  and  cQ-  0,  and  for  y  <  0,  cQ  -  1  and 

Cp-  0,  while  cp  -  0.  As  reactants  are  entrained 

into  the  mixing  cores  of  the  layer,  they  diffuse 
across  the  original  Interface  and  chemical 
reaction  proceeds.  The  rollup  and  pairing 
increases  the  original  length  of  the  Interface  by 
many  folds  and  allow  the  entrained  fluid  to 
diffuse  along  a  larger  boundary  (Ghonlem  et  al. 
[39]).  During  this  process,  if  the  Lagrangian 
elements  utilized  to  represent  the  interaction 
between  chemically  reacting  species  are  brought 
close  enough  so  that  the  distance  between  two 
neighboring  elements  is  smaller  than  the 
characteristic  diffusion  length,  they  react  at  the 
rate  defined  by  Eq.  (17). 

In  Figs.  10,  11  and  12,  we  present  the 
velocity,  location  and  the  strength  of  the 
elements  in  terms  of  product  concentration  for  the 
reacting  mixing  layer  with  constant  rate  chemical 
kinetics  and  temperature-dependent  reaction  rate, 
respectively.  The  amount  of  the  products  formed 
due  to  chemical  reaction  is  represented  by  the 
diameter  of  the  circles  in  the  figures,  l.e. 
larger  circles  indicate  more  products.  In  both 
cases,  Re  »  10000,  Pe  -  4000,  and  U2/U1  •  1/3 
while  Le-1.  In  the  constant  rate  kinetics  case, 
the  value  of  the  Damkohler  number  Da  ■  1  and  In 
the  temperature-dependent  kinetics  Da-  200,  Ta  ■ 
10  and  Q  -  5.  Note  that  in  both  cases  the  value 
of  the  non-dimensional  kinetic  parameters  are  low 
enough  so  that  the  effects  of  heat  release  on  the 
fluid  dynamics  can  be  negligible.  The  stiffness 
of  Eq.  (19)  for  large  values  of  the  Damkohler 
number  Imposes  a  restriction  on  the  time  step  of 
Integration.  In  these  calculations,  we  found  that 
at  -  0.1  is  sufficiently  small  to  accurately 
integrate  the  slow  chemistry . 

A  comparison  between  the  two  figures  reveal 
that  under  Isothermal  conditions,  the  products  are 
formed  as  mixing  occurs  Just  downstream  of  the 
splitter  plate,  while  in  the  temperature-dependent 
kinetics  calculations,  there  is  an  ignition  delay 
before  the  reactant  reach  a  temperature  high 
enough  to  allow  any  significant  chemical  reaction 
to  occur.  Once  the  reaction  begins,  the  mechanism 


Figure  10.  Scalar's  velocity  field  at  Re  *  10000, 
U2/U1  -  1/3- 
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Figure  11.  Product  concentration  field,  Re  ■  10000, 
Pe  •  4000,  U2/U1  »  1/3,  isothermal 
reacting  layer. 


Figure  12.  Product  concentration  field,  variable 
temperature  reacting  layer. 
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of  product  formation  and  chemical  reaction  in  both 
cases  are  asymptotically  the  same.  Increasing  the 
Damkohler  number  to  Da  -  *400  results  in  a  shorter 
ignition  delay,  and  preheating  the  reactants  by 
increasing  the  temperature  at  the  inlet  to  T1  - 
Q/2  while  Da  -  200,  eliminates  the  ignition  delay 
as  indicated  in  Figs.  13  and  14,  respectively. 

In  order  to  examine  the  effects  of  chemical 
reaction  on  the  transport  of  species,  the 
concentration  statistics  in  the  temperature- 
independent  reaction  case  are  presented  in  Figs. 
15  and  16.  These  figures  correspond  to  the 
ensemble  mean  and  fluctuation  in  the  bottom-stream 
species  concentration  in  a  reacting  mixing  layer 
with  Da  •  1,  U2/U1  -  1/2,  Re  -  10000,  and  Pe  - 
4000.  A  comparison  between  figures  15  and  8,  and 
between  figures  16  and  9  indicates  that  near  the 
free  stream,  the  chemistry  affects  the  statistical 
behavior  of  the  species.  Near  the  reaction  zone, 
however,  the  mean  and  the  rms  values  of  the 
concentration  are  lower  under  reacting  conditions, 
while  the  second  hump  near  the  high  speed  stream 
side  of  the  rms  profile  in  the  non-reacting  layer 
is  eliminated  in  the  reacting  flow  due  to  the 
local  consumption  of  the  species  by  chemical 
reaction.  The  same  behavior  was  also  observed  in 
the  experiments  of  Masutani  and  Bowman  [30]  in  a 
reacting  mixing  layer  under  Isothermal  conditions. 
Their  results,  however,  can  not  be  compared 
quantitatively  with  the  present  calculations  since 
the  values  of  the  chemical  parameters  employed  in 
the  numerical  simulation  are  substantially  lower 
than  those  of  the  experiment. 
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Figure  13-  Product  concentration  field,  variable 
temperature  reacting  layer. 
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Figure  14.  Product  concentration  field,  variable 
temperature  reacting  layer. 


Figure  15-  Normalized  mean  concentration  profile 
as  a  function  of  the  cross-stream 
coord inate . 


Figure  16.  Normalized  rms  concentration  profile 
as  a  function  of  the  cross-stream 
coordinate . 


II I. 3  EFFECT  OF  HARMONIC  FORCING 

The  dynamic  effect  of  oscillating  the 
upstream  side  of  the  layer  was  studied 
experimentally  by  several  authors,  e.g.  Oster  and 
Wygnanski  [35]  and  Roberts  and  Roshko  [21]  and 
numerically  by  Ghonlem  and  Ng  [22].  Their  results 
indicate  that  in  the  forced  case,  eddy 
Interactions  follow  four  stages.  In  the  first 
stage,  the  layer  rolls  up  at  the  harmonic  of  the 
forcing  frequency  closest  to  the  most  amplified 
mode.  In  the  second  stage,  a  process  of 
accelerated  pairings  yields  a  large  eddy  which  is 
in  tune  with  the  forcing  frequency.  This  large 
resonant  eddy  appears  earlier  than  it  would  appear 
in  the  case  of  an  unforced  layer.  In  the  third 
stage,  pairing  among  resonant  eddies,  which 
represents  a  neutrally  stable  mode,  is  disabled 
and  the  growth  of  the  vortlcity  layer  is  impaired 
for  several  eddies  downstream.  In  the  fourth 
stage,  the  effect  of  forcing  diminishes  and 
pseudo-random  pairing  is  resumed.  Moreover, 
velocity  statistics  is  affected  by  forcing,  and 
the  sign  of  momentum  transfer  across  the  layer  is 
reversed  following  pairing.  Entrainment  of 
passive  particles  was  found  to  be  commensurate 
with  the  development  of  the  vorticlty  layer. 


In  the  recent  experiment  by  Roberts  and 
Roshko  [21],  it  has  been  observed  that  periodic 
forcing  has  a  direct  influence  on  the  outcome  of 
chemical  reaction  across  a  turbulent  shear  layer. 
The  results  of  this  experiment  indicate  that  when 
harmonic  forcing  is  applied,  the  mixing  rate:  (1) 
is  increased  in  the  Initial  stages  where  the 
resonant  eddy  is  forming:  (2)  is  decreased  in  the 
intermediate  stage  which  corresponds  to  the 
resonant  or  "frequency-locked"  region:  and,  (3)  is 
the  same  as  that  of  the  unforced  layer  further 
downstream.  In  order  to  characterize  these  three 
regions,  the  Wygnanskl-Oster  parameter  X  -  aU  fl 

2  w 
x/Um  is  utilized,  where  ft  is  the  forcing 

frequency  [35].  Roberts  and  Roshko  [21]  and 
Browand  and  Ho  [36]  show  that  the  three  different 
regions  can  be  classified  according  to  the  local 
value  of  Xu  parameter.  In  region  I ,  X  <  1  ,  the 

growth  rate  is  enhanced.  In  region  II,  Xw  >  1,  the 

frequency-locked  region,  the  growth  rate  is 
inhibited.  In  region  III,  the  growth  rate  relaxes 
to  that  of  the  unforced  layer. 

In  order  to  investigate  this  phenomenon 
computationally,  the  response  of  the  reacting 
shear  layer  to  the  application  of  low  frequency, 
low  amplitude  perturbations  on  the  upstream  side 
of  the  shear  layer  is  computed.  Streamwise 
oscillations  are  applied  on  both  sides  of  the 
layer,  hence  a  pressure  perturbation  is  Imposed 
without  changing  the  vorticity  field.  The 
streamwise  velocities  are  taken  as  U1  ■  1  ♦  a  sin 
(2u!lt),  and  U2«  a  U2,  where  a  is  the  amplitude  of 
forcing. 

The  normalized  distribution  of  the  product 
thickness  along  the  mixing  layer  for  three  cases, 
Q  -  0,  0.5  and  1,  is  shown  in  Fig.  17.  In  these 
calculations,  a  -  0.1,  and  Re  -  U000.  The  figure 
Indicates  that  for  n  •  1,  mixing  is  enhanced  in 

the  initial  part  of  the  layer,  1  S  x  S  2.  The 
resonant,  frequency-locked  region  begins  at  x  •  2 
and  ends  at  value  x  -  3.  In  this  region,  mixing 

is  reduced  and  is  less  than  that  of  unforced 
mixing  layer.  Downstream  of  this  region,  x  t  3. 
mixing  rate  resumes  its  natural  growth  and  reaches 
asymptotically  that  of  the  unforced  layer.  For 
lower  forcing  frequency,  G  •  0.5,  the  same  overall 
behavior  is  observed.  In  this  case,  however,  the 
results  of  numerical  calculations  indicate  that 
the  resonant  frequency-locked  region  is 
approximately  in  the  range  3  S  x  S  4.  A 
comparison  between  the  range  of  the  frequency 
locked  region  calculated  here  with  that  estimated 
by  Browand  and  Ho  [36]  is  shown  on  Table  I. 
Considering  the  fact  that  our  simulations  Ignore 
the  effect  of  small  scale  three-dimensional 
turbulence  motion,  and  considering  the  non- 
universality  of  the  Browand  and  Ho's  curve  due  to 
its  Independence  to  experimental  conditions  and 
other  important  non-dimensional lzed  parameters, 
this  agreement  is  encouraging. 


TABLE  I 


Figure  17.  Variation  of  the  product  thickness 
versus  the  downstream  distance. 


Ill  .1)  EFFECT  OF  STRAIN  RATE 

It  has  been  shown  experimentally  by  Tsujl 
[37],  numerically  by  Liew  et  al.  [38],  and 
analytically  by  Peters  [39],  that  the  strain  rate 
has  a  major  influence  on  the  flame  structure, 
particularly  In  non-premixed  systems.  In  the 
counter-flow  diffusion  flame  experiment  of  Tsujl 
[37],  it  was  observed  that  increasing  the 
magnitude  of  stretch  near  the  flame  surface 
results  in  an  increase  of  the  flow  of  reactants 
into  the  reaction  zone.  As  a  result,  the  chemical 
reaction  is  not  able  to  keep  pace  with  the  supply 
of  reactants,  and  the  reaction  rate  is  reduced 
until  local  flame  quenching  occurs.  The  analysis 
of  Peters  [39],  which  is  based  on  the  method  of 
matched  asymptotic  expansion  at  large  activation 
energy,  shows  that  the  mechanism  of  flame 
extinction  can  be  addressed  by  examining  the  local 
value  of  the  rate  of  scalar  dissipation.  This 
parameter  is  viewed  by  Peters  [39]  as  the  Inverse 
of  the  diffusion  time  scale.  If  the  local  value 
of  dissipation  is  Increased  beyond  a  critical 
limit,  the  heat  conducted  away  from  of  the 
diffusion  flame  can  not  be  balanced  by  the  heat 
produced  by  the  chemical  reaction.  As  a  result, 
the  maximum  value  of  the  temperature  decreases, 
and  the  reaction  eventually  ceases. 

By  increasing  the  number  of  scalar  elements 
to  38  in  each  stream  while  decreasing  the 
computational  domain  to  Xmax  •  4,  and  by 
preheating  the  incoming  reactants  to  Ti  -  Q/2  to 
start  the  chemical  reaction  immediately  downstream 
the  splitter  plate,  we  were  able  to  observe  this 
phenomenon.  Figures  18  and  19  show  the 


instantaneous  velocity  and  temperature  rise,  T-Tl, 
of  the  scalar  elements  at  times  of  t  ■  19  and  t  « 
21,  respectively.  In  this  case,  the  Damkohler 
number,  the  normalized  enthalpy  of  reaction,  the 
activation  energy  and  the  velocity  ratio  at  the 
inlet  are  50,  8,  20  and  1/3,  respectively.  The 
cross  stream  direction  is  enlarged  by  a  factor  of 
2  for  the  purpose  of  clarity. 

The  figures  show  that  the  number  of  scalar 
elements  near  the  braid,  which  is  the  thin  link 
between  two  neighboring  cores,  is  only  a  small 
portion  of  the  total  number  of  elements  within  the 
computational  domain,  which  reached  more  than 
5100.  This  indicates  an  Instantaneous  quenching 
at  the  stagnation  points  of  the  layer.  Moreover, 
the  temperature  and  product  concentration  in  the 
reaction  zone  reach  a  maximum  at  the  core  of  the 
eddies  where  the  vorticity  concentration  is  high, 
while  they  reach  a  minimum  at  the  stagnation  point 
within  the  braid  between  the  neighboring  cores 
where  the  strain  and  the  scalar  gradients  reach 
their  maximum  values.  This  is  consistent  with  the 
results  of  the  pseudo-spectral  calculations  of 
Givi  et  al.  [15],  and  with  the  experimental 
observations  of  Tsuji  D7 ]  who  showed  that  the 
local  extinction  of  diffusion  flames  occurs  mainly 
at  the  regions  of  high  dissipation  rate.  At  these 
regions,  the  temperature  tends  to  decrease,  and  if 
it  goes  below  a  critical  characteristic  value,  the 
flame  locally  extinguishes. 

Quantitative  analysis  of  the  effects  of 
stretch  on  the  chemical  reaction  is  rather 
difficult  in  the  context  of  present  algorithm. 
This  is  due  to  the  fact  that  there  are  very  few 
scalar  elements  near  the  regions  of  high  strain, 
and  as  shown  by  Ghoniem  et  al.  [3*0,  most  of  the 
elements  tend  to  be  concentrated  near  the  regions 
with  low  dissipation.  Implementation  of  a 
numerical  scheme  based  on  the  transport  of  the 
scalar  gradients,  as  in  Ghoniem  et  al.  [3*0  can 
Improve  the  accuracy  of  the  analysis 
substantially,  particularly  those  associated  with 
the  effects  of  stretch.  In  this  method,  the 
elements  are  concentrated  near  the  regions  of 
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Figure  19-  Temperature  field  for  reacting  mixing 
layer. 


large  gradients,  or  high  dissipation,  and  hence  a 
smaller  total  number  of  elements  have  to  be 
considered.  The  Implementation  of  this  method  for 
the  numerical  simulation  of  unpremixed  reacting 
flows  is  presently  underway  to  study  the  effect  of 
strain  rate  more  accurately. 


IV  CONCLUSIONS 

In  this  work,  a  numerical  scheme  based  on  the 
transport  of  computational  elements  carrying 
vorticity  and  scalar  quantities  has  been  developed 
to  simulate  a  reacting  planar,  two-stream  mixing 
layer  with  unmixed  reactants.  The  scheme  solves 
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Figure  18.  Temperature  field  for  reacting  mixing 
layer. 


the  transport  equations  at  high  Reynolds  and 
Peclet  numbers  without  using  models  for  turbulence 
closure.  A  Lagranglan  stochastic  model  is  used  to 
implement  the  chemical  reactions  for  both  constant 
rate  kinetics  and  variable  temperature  Arrhenius 
reactions. 

In  the  non-reacting  flow  simulations,  the 
calculated  statistics  of  the  mixing  of  a  conserved 
scalar  are  in  good  agreement  with  experimental 
data.  In  particular,  the  numerical  results  show 
the  presence  or  two  maxima  in  the  fluctuation 
profile.  In  the  constant  rate  reacting  flow 
simulation,  the  effect  of  chemistry  is  to  smooth 
out  this  curve  and  produce  a  single  maximum,  which 
agrees  with  the  experimental  observations. 
Harmonic  forcing  enhances  the  mixing  within  the 
accelerated  growth  zone  of  the  vorticity  layer, 
while  it  impairs  the  entrainment  of  the  unmixed 
fluid  into  the  cores  in  the  resonating  region.  As 
a  result,  the  numerical  simulation  Indicates  a 
decrease  in  the  rate  of  product  formation  in  the 
frequency-locked  region,  similar  to  previous 
experimental  findings. 

In  the  Arrhenius,  temperature-  dependent 
kinetics,  the  mechanism  of  ignition  delay  and  the 
effects  of  reactants  preheating  on  the  decrease  of 
the  duration  of  this  delay  is  observed.  Also,  the 


non-equilibrium  coupling  between  the  scalar 
dissipation  rate  and  the  flame  structure  is 
revealed  as  quenching  frequently  appears  within 
the  braids.  To  describe  this  phenomenon  more 
accurately,  work  is  underway  to  construct  a  higher 
order  scheme  which  can  provide  better  resolution 
at  the  regions  of  strong  strain  rates. 
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The  paper  on  "Three  dimensional  vortex  simulation  with  application  to 
an  axisymmetric  shear  layer"  describes  the  three  dimensional  vortex  element 
method  and  its  application  to  the  evolution  of  the  azimuthal  instability  on 
a  vortex  ring  and  the  initial  stages  of  development  of  a  turbulent  jet. 
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ABSTRACT 

A  three  dimensional  vortex  element  method  is 
developed  for  the  numerical  simulation  of 
incompressible  flow  at  high  Reynolds  number.  The 
method  utilizes  vortex  vector  elements  with  finite 
point-symmetric  cores  to  discretize  the  vortlclty 
field.  The  transport  of  these  elements  is  done  in 
Lagrangian  coordinate  by  computing  the  velocity 
field  as  a  summation  over  the  Individual 
contributions  of  the  elements. 

The  method  is  used  to  compute  the  self- 
induced  velocity  of  a  vortex  ring  and  the 
stability  of  a  vortex  ring  with  finite  core. 
Results  show  that  vortex  rings  become  unstable  to 
a  particular  azimuthal  perturbation  that  depends 
on  the  core/radius  ratio.  The  mode  frequency  and 
shape  of  the  unstable  state  are  in  excellent 
agreement  with  analytical  and  experimental 
results.  The  method  is  applied  to  study  the 
rollup  of  an  axisymmetrlc  shear  layer  and  the 
generation  of  large  scale  vortex  ring  structures. 


I.  INTRODUCTION 

At  high  Reynolds  numbers,  vortlclty  occupies 
a  small  subset  of  the  volume  of  the  flow  field. 
This  is  exemplified  by  boundary  layers,  shear 
layers,  waxes.  Jets,  separation  and  recirculation 
zones,  etc.  These  vorticlty  distributions  are 
unstable  to  natural  perturbations.  At  small 
amplitudes,  perturbations  grow  exponentially  in 
time,  however,  they  have  a  limited  effect  on  the 
flow.  The  growth  of  these  perturbations  into  the 
non-linear  stages  is,  howeve-,  accompanied  with 
severe  distortions  of  the  shape  of  the  vorticlty 
field  and  strong  changes  in  the  local 
concentration  of  vortlclty.  Examples  for  these 
changes  is  the  formation  of  large  scale  structures 
in  shear  layers  and  recirculation  zones. 
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In  highly  three-dimensional  flows,  several 
forms  of  instability  may  arise  simultaneously. 
The  evolution  of  spanwlse  waves  on  the  large  scale 
eddies  and  the  development  of  azimuthal 
instability  along  the  axis  of  vortex  rings  have 
been  observed  experimentally  as  evidence  of 
multiple  forms  of  instability.  In  this  case,  the 
distortion  of  the  vortlclty  field  occurs  faster 
and  the  non-linearity  is  compounded  by  the 
interaction  between  different  instability  modes. 
Moreover,  the  problem  is  governed  by  several 
length  and  time  scales,  and  multiple  states  can  be 
expected  depending  on  which  mode  grows  faster  (for 
photographic  record  of  the  development  of 
vorticlty  fields,  see  Van  Dyke  [1]  and  Lugt  [2]). 

It  has  been  reported  experimentally,  and 
observed  in  numerical  studies,  that  these  changes 
in  the  vortlclty  field  may  not  incur  strong 
variations  in  the  swan  flow  rield.  This  is 
expected  since  the  velocity  Is  an  Integral  mean  of 
the  vortlclty  field.  However,  they  arfect  the 
fluctuations  strongly  and  to  the  level  where  the 
order  of  magnitude  of  the  fluctuation  may  change. 
This  is  extremely  Important  in  mixing  and  heat 
release  in  chemically  reacting  flows  since  the 
rate  of  mixing,  and  thus  chemical  reaction,  is  a 
strong  function  of  the  fluctuation  and  depends 
weakly  on  the  mean  field.  It  has  been  confirmed 
that  by  changing  the  vortlclty  field  of  a  shear 
layer  through  imposing  certain  perturbations  on 
the  flow,  the  rate  of  chemical  reaction  can  be 
enhanced  or  slowed  and  that  turbulent  shear 
stresses  can  reverse  sign  during  the  same  process 
(for  a  review  and  some  recent  results,  see  Ho  snd 
Huerre  [3],  Robert  and  Roshko  [k]  and  Ghoniem  and 
Ng  [53). 

To  capture  these  changes,  numerical 
simulation  of  the  unaveraged  non-linear  aquations 
of  motion  has  been  utilized.  For  the  success  of 
these  simulations,  care  must  be  exercised  in 
resolving  small  variations  since  they  ultlmstely 
grow  to  produce  the  finite  amplitude  changes,  and 
hence  numerical  diffusion  should  be  minimized. 
Moreover,  schemes  must  adapt  to  the  strong 
distortion  in  the  flow  field  without  developing 
numerical  Instabilities.  Thus,  Lagrangian  schemes 
seem  like  natural  candidates.  A  grid-free  class 
of  Lagrangian  schemes,  vortex  methods,  is  utilized 
in  this  work  to  study  the  evolution  of  three 
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dimensional  flow  fields  dominated  By  high 
concentration  of  vorticity  at  n.ign  Reynolds 

numbers. 

In  vortex  simulation,  vorticity  distribution 
Is  represented  by  a  finite  number  of  localized 
vortex  elements,  or  vortex  disks  In  two  dimensions 
(Blobs  according  to  Chorin  [6]),  and  vortex  balls 
In  tnree  dimensions  (vortons  according  to  Saffman 
and  Merlon  [7],  or  vortex  arrows  according  to 
Leonard  [8]'  which  move  in  an  inviscld  field  with 
the  local  particle  velocity.  A  particular  flow 
conf lguratlon  can  be  fully  described  when  the 
appropriate  boundary  conditions  are  Imposed  on  the 
velocity  field  by  adding  an  extra  irrotatlonal 
field.  Two-dimensional  vortex  simulations  have 
been  useful  In  providing  an  accurate  description 
of  the  large-scale  structure  of  turbulence  in 
shear  flows  (Ashurst  [9],  Ghonlem  and  Sethian  [10] 
and  Ghonlem  and  Ng  [5]).  However,  they  cannot  be 
used  to  describe  phenomena  In  which  streamwlse 
vorticity,  or  variation  along  vortex  lines,  plays 
an  Important  role.  Moreover,  they  lack  the 
ability  to  capture  small-scale  turbulence 
structures  which  arise  due  to  vortex  stretching 
and  tilting  with  reaped  to  the  main  flow  plane. 

In  this  work,  a  three  dimensional  vortex 
element  method  is  developed  for  the  numerical 
simulation  of  flow  field  with  high  concentration 
of  vorticity  at  high  Reynolds  number.  The  scheme 
utilizes  vortex  vector  elements  with  finite  point- 
symmetric  cores  to  discretize  the  vorticity  field 
and  follows  their  motion  in  Lagrangian 
coordinates.  The  vortex  vector  elements  change 
their  vorticity  according  to  the  local  stretch, 
while  their  direction  is  determined  By  the  tilting 
of  the  vortex  lines.  The  rotational  velocity 
field  is  computed  by  summing  over  the  field  of 
each  individual  element,  which  is  evaluated  from 
tne  desingularized  Biot-Savart  law.  The  potential 
velocity  added  to  satisfy  the  boundary  conditions 
is  computed  by  using  the  appropriate  image  system 
of  the  vortices.  For  recent  reviews  of  vortex 
calculations  In  three  dimensions,  see  Leonard 
[8,li]  and  Saffman  and  Baker  [12], 

To  check  the  accuracy  of  the  vortex  method, 
we  use  test  problems  and  make  comparison  with 
experimental  and  analytical  results.  The 
discretization  algorithm  is  applied  to  compute  the 
self-induced  velocity  of  a  vortex  ring  and  the 
results  are  compared  with  the  Saffman’s  analytic 
solution  [13].  The  stability  results  of  a  vortex 
ring  with  a  finite  non-deformable  core  and  a 
vortex  torus  with  a  deformable  core  are  compared 
with  the  analytical  solutions  of  Hldnall  and 
Sullivan  [ik]  and  Wldnall  et  al.  [15]. 
Preliminary  results  fo*-  the  rollup  of  a  three 
dimensional  shear  layer  subject  to  an  axl- 
symmetrlc  perturbation  are  compared  with  the 
experimental  results  of  Vandsburger  et  al.  [16] 
and  Roquemore  et  al.  [17]. 


II.  FORMULATION  AND  NUMERICAL  SCHEME 

In  this  section,  the  construction  of  a  three- 
dimensional  scheme  for  tracking  the  evolution  of  a 
vorticity  structure  in  an  arbitrary  domain  Is 
described.  The  accuracy  of  the  scheme  is  checked 
against  theoretical  results  regarding  the  motion 


and  stability  of  a  vortex  ring  and  a  vortex  torus. 
The  scheme  is  constructed  as  follows:  The 
vorticity  field  is  first  discretized  into  a  finite 
number  of  small  straight  line  vortex  vector 
elements  and  then  followed  in  a  Lagrangian  frame 
of  reference.  Each  element  has  a  finite  core  of 
vorticity  which  is  point-sywnetrleal  around  its 
center,  and  hence  the  nomenclature  "vortex  ball". 
The  velocity  produced  by  a  distribution  of  vortex 
vector  elements,  or  vortex  balls,  is  obtained  by 
the  desingularized  Biot-Savart  law,  which  amounts 
to  summing  the  velocity  produced  by  individual 
vector  elements.  The  procedure  for  a  consistent 
discretization  and  the  evaluation  of  the  Biot- 
Savart  law  is  explained  in  Section  II. 2.  Its 
numerical  accuracy  and  convergence  under  steady 
state  conditions  is  shown  in  Section  II. 3.  The 
comparison  between  the  numerical  and  analytical 
results  for  the  stability  a  thin  vortex  ring  and  a 
vortex  torus,  another  test  for  the  accuracy  of  the 
method  under  unsteady  state,  is  discussed  in 
Sections  II. k  and  II. 5. 

The  potential  velocity  field  added  to  satisfy 
a  particular  set  of  conditions  on  the  boundaries 
is  determined  by  solving  the  Laplace  equation 
nume-lcally  subject  to  the  appropriate  boundary 
conditions.  When  the  boundary  conditions  match 
those  of  a  standard  potential  solver,  l.e. 
Dlrichlet  or  Neumann  conditions,  that  particular 
routine  can  be  used  to  evaluate  the  potential 
velocity.  In  cases  when  the  boundary  conditions 
are  neither  Dlrichlet  nor  Neumann  type,  one  faces 
difficulty  in  satisfying  continuity  along  the 
boundaries,  and  a  special  algorithm  must  be 
constructed  to  handle  this  difficulty.  This  is 
discussed  in  Section  III, 

1 1.1.  EQUATIONS  OF  MOTION 

The  motion  of  an  Incompressible  inviscld  flow 
is  governed  by  the  Euler  equations: 


7  •  u  -  0 

(1) 

*  u  ■  Vu  -  -  Vp 

(2) 

expressing  the  conservation  of  mass  and  momentum, 
respectively.  In  these  equations,  u  •  (u,v,w)  is 
the  velocity,  t  is  time,  ?  •  (3/3x,3/3y,3/3z)  is 
the  gradient  operator,  while  1  •  (x,y,z),  and  p  is 
pressure.  Quantities  are  normalized  with  respect 
to  the  appropriate  combination  of  a  characteristic 
length  scale,  velocity  scale  and  density.  In 
vortex  simulation,  the  equations  are  recast  in 
terms  of  the  vorticity  «: 

•  •  ?  x  u  (3) 

The  vortex  transport  equations  are  obtained 
by  taking  the  curl  of  Eq.  (2).  Using  Eq.  (1)  and 
using  the  fact  that  V.w  •  V.YXu  -  0,  l.e.  the 
vorticity  rorms  a  aolenoldal  vector  field,  we  get: 

♦  u  •  fa  -  *  •  fu  (A) 

If  the  vorticity  field  Is  known,  the  velocity  can 
be  evaluated  by  integrating  Eqs.  (1)  and  (3),  as 
shown  Below,  while  Eq.  (A)  is  used  to  transport 
the  vorticity  in  the  form  of  a  number  of  discrete 
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elements.  This  Is.  In  essence,  the  backbone  of 
vortex  simulation  algo-1  thins. 

Based  on  the  Helmholtz  decomposition  of  a 
vector  field,  the  velocity  can  be  spilt  Into  a 
solenoldal  and  an  irrotatlonal  component, 

u  *  u  *  u  (5) 

U)  P 

where  u  Is  the  velocity  Induced  by  the  vortlclty 
u 

field  In  an  unbounded  space,  and  up  Is  the 

potential  component  added  to  satisfy  the  potential 
no-flow  condition  along  the  boundary  of  the 
domain.  Each  component  will  be  evaluated 
separately,  satisfying  the  appropriate  boundary 
conditions,  and  then  added  to  obtain  the  total 
velocity. 

To  evaluate  the  velocity  field  induced  by  a 
given  vorticity  distribution  ■  In  an  unbounded 
space,  we  assume  the  existence  of  a  vector  stream 
function  f  such  that 


Furthermore,  the  vortlclty  vector  can  be 
written  as  ■  -  u  I.  where  1  Is  the  direction  of 
the  local  vortex  line,  or  material  lines,  while  dx 
•  dA.dA  where  A  Is  the  cross  sectional  area  normal 
to  the  direction  A.  The  circulation  of  a  vortex 
line,  r,  which  is  conserved  along  a  particle  path 
according  to  Kelvin  theorem,  Is  expressed  in  terms 
of  the  vortlclty  »  as  r  •  (/  m.dA. 

Since  u  Is  an  Irrotatlonal  component,  then 
P 

u  •  9e.  where  p  Is  a  velocity  potential  governed 

P 

by: 

V2*  .  0  O' ) 


while  the  total  velocity,  u.  Is  subject  to  a 
potential  boundary  condition  at  the  boundaries  of 

the  domain,  l.e.(u  ♦  ¥p).n  Is  given  on  3D,  where  n 
u 

Is  the  local  normal  to  3D,  which  denotes  the 
boundary  of  the  domain. 


u  «  V  x  p  (6) 

w 

By  construction,  u  satisfies  the  continuity 

w 

equation  since  V.¥xp  -  0  Identically. 

Substituting  In  Eq.  (3)  and  assuming  tnat  V.p  *  0, 
one  obtains: 

V2  «  -  -  M  (7) 


The  solution  of  this  Poisson  equation  In  three 
dimensions  Is  given  by: 

*(x)  -  !  C(x  -  x‘ )  mix')  dx'  (8) 

where  G(x)  -  1/<&wr)  Is  the  Green  function,  and  r 
•  ) x l .  As  shown  by  Batchelor  [18],  the  above 
expression  for  p  Is  solenoldal,  as  previously 
assumed,  if  the  boundaries  of  the  domain  extend  to 
Infinity.  This  Is  essentially  the  condition  needed 
to  evaluate  u  . 

The  solenoldal  velocity  component,  u^.can  be 

evaluated  by  substitution  in  Eq.  (6)  which  yields 
the  well-known  Biot-Savart  law: 

u  •  /  K'x-x')  x  mix')  dx'  (9) 


K(x) 


(10) 


where  1'  is  the  position  of  the  volume  element 
dx’ . 

The  implications  of  the  equations  of  motion, 
Eqs.  (3)  and  ( u ) ,  regarding  the  evolution  of  the 
vortlclty  field  can  be  summarized  In  the  following 
Important  dynamical  statements,  given  here  without 
proof  while  used  later  in  the  construction  or  the 
numerical  algorithm  (for  details,  see  Batchelor 
[18]): 

(1)  Kelvin  theorem:  The  circulation  around  a 
closed  material  loop,  defined  as  r  •  Af  m.dA  where 

A  Is  the  surface  area  within  the  loop,  Is 

conserved  as  the  loop  Is  deformed; 

(2)  Helmholtz  theorem:  Vortex  lines, 

parallel  everywhere  to  the  local  vortlclty  vector, 
move  as  material  lines. 


II. 2.  EVALUATION  OF  THE  ROTATIONAL  FIELD 

Analytical  evaluation  of  the  Biot-Savart 
integral  in  Eq.  (9)  Is  restricted  to  simple 
vortlclty  distributions.  such  as  rectllnear 
vortices  and  circular  vortex  rings  of  concentrated 
vortlclty.  Therefore,  the  Integration  must  be 
performed  numerically  for  an  arbitrary  vortlclty 
distribution.  For  that  purpose,  the  continuous 
vortlclty  field  is  discretized  Into  a  number  of 
vortex  vector  elements,  each  with  an  assigned 
vortlclty  Mj.  The  magnitude  of  the  vortlclty 

associated  with  each  element  Is  distributed  over  a 
small  spherical  volume  around  Its  center  according 
to  a  core  function  f  with  a  characteristic  core 
radius  4.  The  vortlclty  field  Is  hence  expressed 
as: 


»(x,0) 


N 
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•t(0)  h3  f4(x  -  x4) 


(12) 


where  N  Is  the  total  number  of  vortex  vector 
elements,  and  Xj  Is  the  center  of  the  vortex 

element,  while  h  Is  the  Initial  distance  between 
the  centers  of  neighboring  elements.  The  accuracy 
of  this  discretization  Is  discussed  in  Beale  and 
Majda  [19,20].  Note  that  If  f{-  A(x-Xj),  where 

S(.)  Is  the  Dirac  delta  function,  Eq.  (12) 
represents  a  distribution  of  singular  vortex 
lines,  or  vortons  [7].  However,  in  this 
representation,  vortex  balls,  while  equivalent  to 
vortex  disks  In  two  dimensions,  are  used  and  f.- 

i  6 

I/43  f(r/4),  while  6  is  finite.  The  distribution 
of  the  magnitude  or  the  vorticity  associated  with 
each  element  Is  point  symmetrical  around  Its 
center  jj,  while  Its  direction  everytdiere  Is  l- 

mj/wj,  and  w  -  |m|.  «  Is  the  core  radius  of  the 

element  where  most  of  Its  vorticity  Is 
concentrated,  f  >  0  for  r  <  6,  and  f  vanishes 

rapidly  for  r  >  6. 

A  simple  Intuitively  appealing  choice  for  a 
core  function  could  be  the  Hill  spherical  vortex 
for  which  u  •  A  p  for  p  <  4  and  u  •  0  for  p  >  4, 
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wnere  p2  -  x2  ♦  y2,  o  being  the  radial  direction 

In  a  cylindrical  coordinate  system.  This, 
however,  Is  a  poor  choice  for  the  core  function 
since  the  latter  should  have  a  maximum  at  the 
center  and  decay  further  away.  Moreover,  u  la  not 
a  function  of  r.  A  better  choice  may  be  the  three 
dimensional  analogue  of  the  Rankine  vortex,  l.e.  « 
•  w-  for  r  S  a  and  y  •  0  for  r  J  i. 

For  accurate  discretization  of  the  continuous 
vortlclty  distribution,  6  should  be  chosen  larger 
than  h,  where  h  Is  the  Initial  separation  between 

the  elements  and  h3  Is  the  volume  element  used  to 
construct  the  vortex  balls.  This  will  ensure  that 
the  core  functions  f  associated  with  neighboring 
elements  are  highly  overlapping.  Tne  Introduction 
of  a  similar  discretization  procedure  has  been 
widely  used  in  two  dimensional  simulations  to 
construct  stable  and  accurate  vortex  algorithms 
[5,6,8,9,19],  Moreover,  h  may  take  on  different 

values  In  the  three  principle  directions,  and  h3 
is  replaced  by  oV,  where  av  -  h^  n^  h^.  Ir.  this 
case,  the  vortlclty  associated  with  an  element  is 
•  1/  flV  /  u  dx,  the  integration  Is  performed 
over  aV. 

From  Helmholtz  theorem,  the  vortlclty 
associated  with  a  material  element  alj  changes  as 

It  stretches,  «  '  t)  •  l*»,  ( 0)/att<  0) )  at^Ctl,  while 

at  *  |al|.  Moreover,  according  to  Kelvin  theorem, 
the  circulation  of  the  vortex  vector  element 
remains  constant  as  It  moves  along  particle  paths, 
while  due  to  lncompressblllty.  a v  Is  constant. 
Thus,  Eq.  (11)  can  be  written  as: 

N 

m(x,t)  -  l  Ti  alj(t)  rfi( *~xi >  03) 


Thus,  using  a  first-order  time  Integration! 

g^t-at)  *  gjtt)  ♦  Uj  at  07) 

and 

at^trat)  •  atjU)  ♦  al^t)  •  fu^  at  08) 

The  velocity  gradient  can  be  evaluated 
analytically  by  differentiating  Eq.  0*0.  •* 

proposed  by  Anderson  and  Greengard  [22].  However, 
since  the  vortlclty  moves  along  particle  paths, 
the  material  line  coinciding  with  the  vortex 
vector  alj  can  be  approximated  by  lta  terminal 

points  *it-  ♦  alj/2  and  *2^  »j  -  alj/2.  and 

the  center  of  the  vortex  vector  element  Is 
approximated  by  g  -  (gi  ♦  l2)/2.  In  this  scheme, 
a  vortex  vector  element  Is  described  by  (r,  gi , 
g2)  and  both  terminal  points  are  updated  each  time 
step.  A  similar  construction  was  used  by  Chorin 
[2 3, 2 A, 25]  to  study  boundary  layer  stability,  the 
evolution  of  a  turbulent  vortex  and  the  properties 
of  developed  turbulence.  Since  the  vortlclty 
field  Is  adenoidal,  the  end  of  an  element  is  the 
beginning  of  the  next  element  If  these  elements 
were  neighboring  elements  on  the  same  vortex  line 
at  t  •  0.  Thus,  this  scheme  can  be  used  to  ensure 
the  satisfaction  of  the  condition  V.m  •  0  by 

maintaining  the  connectivity  of  vortex  tubes  no 
matter  how  accurate  is  the  discretization  or  the 
vortlclty  field.  Tne  same  property  Is  utilized  In 
the  filament  scheme  of  Leonard  [26]  (see  also 
Ashurst  and  Melburg  [27]).  A  discussion  of  the 
relationship  between  different  algorithms  is  given 
by  Greengard  [28].  In  our  computations,  a  second- 
order  time  integration  algorithm  is  used  to  move 
the  terminal  points,  (gi.g2).  of  the  vortex  vector 
elements,  e.g.,  for  gi : 

gi*  -  4 


In  this  expression,  at,  Is  the  material  vector 
associated  with  the  vortix  vector  element,  and  g 
is  the  midpoint  of  this  vector,  g  (Xj.O)  •  Xj. 

The  velocity  field  Is  obtained  by 
substituting  Eq.  (9)  Into  Eq.  (10)  and 
integrating: 

Uu  *  -  k 
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where  k(r)  -  A*  /  f(r’)  r'2  dr’,  and  ^  •  |x  - 
gj,  while  xl  is  the  center  of  the  vortex  vector 
alj.  In  this  representation,  each  vortex  vector 
element  is  described  by  (r,  g,  aty  Since  y  and 

at  are  the  position  and  length  of  a  material 
particle  and  a  material  line,  respectively,  their 
variation  with  time  can  be  obtained  from 
(Truesdell  [21]): 


U  *  ■c«(x.t>.t> 


05) 


gij(t»at)  •  *ij(t)  ♦  (Uj  ♦  Uj)/2  at 
where  u*  •  u(g*). 

The  accuracy  of  the  vortlclty  discretization 
depends  on:  the  choice  of  the  core  function  f,  the 
distance  between  the  centers  of  neighboring 
elements  h,  and  the  ratio  between  the  Initial 
separation  between  the  vortlcltles  and  the  core 
radius,  «/h.  In  tne  analysis  of  Beale  and  Majda 
[29],  It  Is  shown  that  a  second  order  scheme  Is 
obtained  If  the  following  third  order  Gaussian 
core  function  Is  used: 

3  -r3 

f(r)  -  {7  •  (20) 

and; 


<(r)  -  1  -  e"r  (21 ) 


dal 

dt 


at  •  ?u 


(16) 


A 


As  the  flow  develops  strong  stretch  along  the 
vortex  lines,  the  effective  value  of  at  exceeds  h 
and  the  amount  of  vortlclty  transported  by  each 


vortex  bell  grows.  To  maintain  a  uniform 
resolution,  if  atj  >  2h,  a  vortex  element  is  split 

Into  two  elements,  each  with  at  •  atj/2  and  r  • 
fj.  Effectively,  this  amounts  to  redistributing 

the  vorticlty  field  among  a  larger  number  of 
elements  to  maintain  the  accuracy  of  the 
calculations.  The  value  of  6  is  kept  constant  in 
our  calculations. 

II. 3.  SELF-INDUCED  VELOCITT  OF  A  SING 

To  investigate  the  effect  of  discretization 
of  the  vorticlty  field  on  the  accuracy  of  the 
calculation  of  the  velocity  field,  the  self- 
induced  velocity  of  a  vortex  ring  with  a  radius  R 
and  a  finite  core  radius  o  is  computed  and 
compared  with  the  analytical  results  for  a  thin 
vortex  ring.  For  this  purpose,  the  ring  is 

discretized  along  its  axis  into  a  number  N  of 
vortex  vector  elements,  where  the  length  of  each 
element  is  h  •  at(0)  -  2»R/N.  Each  element  is 

represented  by  a  computational  vortex  ball  with 
core  radius  t  •  a.  This  is  a  worst-case  analysis, 
since  normally  one  would  use  several  elements  to 
represent  the  core,  and  choose  6  <  o  (as  will  be 
shown  later).  However,  we  start  with  this  case 
for  simplicity  and  computational  convenience. 

To  distinguish  between  the  two 
representations  of  a  vortex  ring;  where  the  vortex 
balls  are  aligned  along  the  ring  axis  forming  a 
tube  of  vorticlty,  is  called  the  thin  tube 
approximation  while  if  several  vortex  balls  are 
used  within  the  cross  section  of  the  ring,  it  is 
called  a  vortex  torus.  The  first  approximation  is 
different  from  the  thin  filament  approximation  of 
Leonard  [26].  In  the  thin  filament  approximation, 
the  3iot-Savart  law  is  modified  to  remove  the 
singularity  at  the  center  of  the  filament  and  the 
core  maintains  its  vorticlty  distribution  aa  the 
filament  is  deformed.  In  the  thin  tube 

approximation,  neighboring  elements  can  move 
freely  with  respect  to  each  other,  and  hence 

changing  the  local  vorticlty  distribution  of  the 
tube. 

In  the  discretization  of  the  vortex  ring 
using  the  thin  tube  approximation.  6  *  8  h,  where 
8  >  1  to  Insure  the  overlapping  between 

neighboring  elements.  Eq.  (ID)  is  used  to 
evaluate  the  self-induced  velocity  V  by  summing 
the  contribution  of  the  elements  around  the  ring, 
excluding  the  effect  of  the  element  on  itself. 

Results  are  compared  with  the  analytical 

expression  of  Saffman  [13]  for  a  thin  vortex  ring, 

o/R  «  1 : 

V*I7R(ln?-C)  (22) 

where  C  •  0.558  for  a  second  order  Gaussian 

distribution  of  vorticlty  within  the  core  and  0  is 
the  effective  radius,  l.e.  it  is  the  standard 
deviation  in  the  Gaussian. 

A  comparison  between  our  computations  of  the 
self-induced  velocity  and  Eq.  (22)  is  shown  in 

Fig.  1  for  different  values  of  N,  V  •  V/tr/AvR). 
Three  comments  should  be  made  here:  (')  since  the 


N 


Figure  1.  The  normalized  self-induced  velocity  of 

a  vorter.  ring  V  -  V/(r/A»R)  vs.  the  number  or 
computational  vortex  balls  used  to  discretize  the 
ring,  N.  The  analytical  results  of  Saffman  Is 
represented  by  the  straight  line.  o/R-0.1  *  o  ; 
s/R*0.2  *  *  ;  o/R-0.3  *  v. 


core  function  of  the  vortex  elements  is  a  third 
order  Gaussian,  Eq.  (18),  and  not  a  second  order 
Gaussian  as  in  Saffman's  calculations,  a  slight 
discrepancy  in  the  self-induced  velocity  is 
expected  (a  comparison  between  the  two 
distributions  is  presented  in  Fig.  2);  (2)  since  6 
>  h,  and  a  strong  overlap  between  the  cores  of 
neighboring  elements  is  ensured,  the  vorticlty  at 
any  point  is  the  contribution  of  many  elements 
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Figure  2.  A  comparison  between  the  vorticlty 
distribution  within  the  core  of  the  thin  ring  used 
in  the  computations  of  Fig.  1,  l.e.  a  third  order 
Gaussian  described  by  Eq.  (18),  and  that  of 
Saffman  ring,  i.e.  a  second  order  Gaussian.  In 
both  cases,  o/R  •  0.25. 


•Ions  the  ring  axis  (the  vortlclty  of  a  single 
element  and  the  total  ring  vortlclty  at  any  point 
Is  presented  in  Fig.  3,  showing  that  while  the 


Figure  3-  The  vortlclty  distribution  of  a  single 
vortex  ball  and  the  vortlclty  distribution  of  a 
vortex  ring,  both  normalized  with  respect  to  their 
maxima  a  •  6  •  0.25  R. 


magnitude  Is  strongly  affected  be  neighboring 
elements,  the  core  size  Is  the  same);  and,  (3)  the 
analytical  expression  applies  for  a  <<  R,  and 
hence  beat  comparison  is  expected  for  a/R  •  0.1, 
while  It  deteriorates  for  thicker  rings. 

As  the  ring  becomes  thinner,  l.e.  larger  R/o, 
more  elements  are  required  to  achieve  an  accurate 
discretization.  This  Is  expected,  since  by 
choosing  6  •  a  and  6  •  Bh,  where  B  Is  a  factor 
larger  than  one,  the  number  of  elements  N  - 
2*R/(8h)  -  0(R/o! ,  which  increases  as  a  decreases. 
Therefore,  for  a  fixed  core  size  a,  the  number  of 
elements  required  to  compute  the  self-induced 
velocity  due  to  curvature  R  grows  as  R  Increases. 

II. u.  STABILITY  OF  A  THIN  VORTEX  RING 

A  more  interesting  problem,  providing  a  test 
for  the  accuracy  of  the  time-dependent 
computations,  is  the  growth  of  small  perturbations 
on  the  vortex  ring.  There  Is  a  rigorous 
asymptotic  linear  theory  for  the  stability  of 
vortex  rings  In  two  forms:  (1)  for  a  ring  with  a 
non-deformable  core,  performed  by  Wldnall  and 
Sullivan  [lit]  ;  and  (.2)  a  more  elaborate  theory 
for  a  ring  with  a  deformable  core  reported  In 
Wldnall  et  al.  [15]).  Wldnall  [30]  and  Wldnall  ana 
Tsai  [31].  We  will  compare  the  results  of  the 
thin  tube  approximation  to  the  first  case,  and 
results  for  the  vortex  torus  to  the  second  case. 

To  study  the  linear  stability  of  thin  vortex 
rings  in  the  thin  tube  approximation,  l.e.  with 
almost  non-deformable  cores,  a  radial  perturbation 
of  amplitude  c/R  ••  0.02  and  a  wavenumber  n  is 
imposed  on  the  axis  of  the  ring.  The  wavenumber 
is  defined  here  as  the  number  of  waves  that  is 
fitted  along  the  entire  length  of  the  ring  axis. 


tnus  the  size  of  the  perturbation  varies  in  the 
azimuthal  direction  as  ap  ■  c  aln  (2«ne). 
Originally,  the  ring  lias  in  the  x-y  plane,  and 
the  streamwlse  is  the  z-dlrectlon  while  p  •  R  for 
all  vortex  balls.  We  start  with  n  •  1  and 
Increase  by  an  Increment  of  one.  The  time  step 
used  Is  at  •  0.10.  Similar  results  were  obtained 
for  p/R  •  0.1,  0.15,  0.2,  0.25.  In  the  following, 
only  the  first  case  Is  discussed  In  detail. 

For  n  <  n  ,  where  n  Is  the  wavenumber  of  a 
n  n 

neutrally  stable  mode  that  neither  rotates  around 
the  ring  axis  nor  grows,  the  waves  rotate,  or 
spin,  around  the  ring  axis  at  a  frequency  Q  that 
depends  on  n.  As  it  rotates  around  the 
unperturbed  axis  of  the  ring,  the  Instantaneous 
center  of  the  ring  draws  an  ellipse  whose  major 
axis  is  in  the  radial  direction,  p,  and  the  minor 
axis  Is  In  the  streamwlse  direction  z.  Thus, 
these  are  bending  waves  that  move  around  the  ring 
axis,  hence  the  name  helical  waves  (if  the  ring  is 
opened  to  form  a  rectilinear  vortex,  the  waves 
will  like  a  corkscrew  spinning  at  frequency  Q). 
The  sense  of  rotations  of  the  waves  la  the  same  as 
that  of  ring  vortlclty.  The  frequency  of  rotation 
a  starts  out  low  at  small  n,  grows  to  a  maximum 
and  then  decreases  again.  The  amplitude  In  the 
radial  p-dlrectlon  and  streamwlse  z-dlrectlon  are 
shown  in  Fig.  «  for  n  «  2,  A,  6,  8,  10  and  12  for 
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Figure  5.  The  amplitude  of  the  perturbation  in  the 
radial  and  streamwlse  directions  for  the  neutrally 
stable  mode  for  the  same  ring  as  in  Fig.  *. 


produced.  Moreover,  no  wave  rotation  is  observed. 
The  wave  amplitudes  are  shown  in  Fig.  6  for  o/R  • 


Tift  (NU  -  12) 

Figure  «.  The  amplitude  of  the  perturbation  in  the 
•■adial  (('direction  and  the  streamwlse  z-direction 
for  a  vortex  ring  with  o/R  -  0.',  computed  using 
the  thin  tube  approximation.  the  wavenumber  of 
the  perturbation  is  n  -  2.  n,  6,  8,  10  and  12, 

arranged  from  the  top  figure.  Both  amplitudes  are 
normalized  with  respect  to  the  initial 
perturbation  In  the  radial  direction,  t/R  0.02, 

and  time  is  normalized  with  respect  to  R?/r.  In 
this  figure,  the  behavior  of  the  modes  n  <  n  Is 
shown. 
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Figure  6.  The  amplitude  of  the  perturbation  in  the 
radial  and  streamwlse  directions  for  the  unstable 
mode  n*  of  the  ring  of  Fig. 


o^R  *  0.1.  Note  that  the  radial  perturbation 
produces  a  streamwlse  perturbation  of  almost  the 
same  magnitude.  All  these  modes  are  characterized 
as  being  linearly  stable  since  their  amplitudes 
remain  bounded. 

At  n  •  nn,  the  waves  neither  grow  nor  rotate. 
For  o/R  •  0.1,  it  n  •  13  the  ring  remains  In  Its 

original  plane  without  bending,  as  depicted  by 
Fig.  5.  The  next  mode,  n*  •  14,  the  waves  grow  in 
the  radial  direction,  and  then  in  the  streamwlse 
direction  so  that  the  total  amplitude  grows 
exponentially  In  time,  l.e.  the  ring  becomes 
linearly  unstable  and  streamwlse  vortleity  is 


As  n  >  n»,  the  ring  is  stabilised  again  and 
the  eigenfunctions  behave  in  a  similar  way  to  n  < 
nn.  However,  the  major  axis  of  the  ellipse  is  now 

in  the  streaawise  direction  and  the  frequency  of 
rotation  increases  Indefinitely.  Moreover,  the 
sense  of  rotation  Is  In  the  opposite  to  that  of 
the  ring  vortleity.  The  wave  amplitudes  in  the  p 
and  z  directions  are  shown  in  Fig.  7  for  n  •  15, 
U  and  19. 

Similar  observations  are  made  for  o/R  *  0.15, 
0.2  and  0.28.  In  all  cases,  the  unstable  mode  n* 
is  a  bifurcation  in  the  eigenfunction  that 
corresponds  to  0  •  0.  In  Fig.  8,  the  results  of 
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Figure  '.  The  amplitude  of  the  perturbation  In  the 
radial  and  streamwize  directions  Tor  the  modes  n  > 
n*  toe  ring  described  in  Figure  M.  The 
wavenumber  of  the  perturbation  is  n  -  15,  17,  and 
’9.  arranged  from  the  top. 


these  computations  are  summarised  in  terms  of  Q, 
the  rotation  frequency  of  the  waves,  v.s.  the  wave 
number  k  •  no/R.  In  this  figure,  Q  Is  normalized 

2 

with  respect  to  r/R  .  In  all  cases,  the  unstable 
mode  k*  •  n*a/R  *  1.25  corresponds  to  a  non- 
rotating  mode,  Q  •  0.  This  is  in  agreement  with 
the  theoretical  results  of  Widnall  and  Sullivan 
[i«].  They  observed,  similar  to  what  we  see  in 
the  nurne-tcal  results,  that  a  mode  becomes 
unstable  when  the  self-induced  rotation  of  the 
waves  balances  the  rotation  induced  by  the  rest  of 
the  ring  and  the  energy  of  the  perturbation  is 
utilized  in  stretching  the  wave  amplitude. 

In  order  to  check  on  the  accuracy  of  the 
computations,  we  varied  the  discretization 
parameter  h  by  using  more  elements  around  the  axis 
of  the  ring,  h  •  2«R/N.  Figure  9  shows  the  growth 
of  the  amplitude  of  the  perturbation  with  time 
using  an  Increasing  number  of  elements  for  e/R  • 


K 


Figure  8.  The  computed  results  for  the  dispersion 
relation  of  a  ring  using  the  thin  tube 
approximation.  s/R»0.i  *  v  ;  o/R-0.15  ♦  •  ; 

e/R»0.2  -  o  i  s/R-,25  -»o.  The  frequency  of 
rotation  of  the  mode  0  is  normalized  with  respect 

to  R^/f  while  k  Is  normalized  with  respect  to  R/o. 
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Figure  9.  The  growth  of  the  unstable  mode  n»  .  7 
for  o/R  •  o.l,  computed  using  N  •  30-1*0,  with 
Increments  aN  •  10. 


0.2  and  n  ■  n»  .  7.  N  •  30  is  the  smallest  number 
necessary  to  satisfy  the  oondltlon  a  t  h,  however, 
we  notice  that  N  •  90  Is  necessary  to  compute  the 
logarithmic  growth  rate  accurately.  It  is  the 

same  number  required  to  compute  V  •  3.1309 
accurately,  as  seen  in  Fig.  2,  This  Is  not 
surprising  since  the  stability  of  the  local  waves 
depends  strongly  on  the  velocity  (or  strain  field) 
induced  by  the  ring  on  the  perturbation.  The 
linear  growth  rate,  •  a(log«)/at  •  0.1625.  The 
analytical  value  obtained  by  Kidnall  and  Sullivan 
['*]  for  the  same  value  of  V  it  .  0.157. 
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Figure  10.  Tne  computed  wavenumber  of  the  most 
unstasle  mode  vs.  tne  normalized  self-induced 

velocity  V,  compared  wltn  the  results  of  the 
linear  asymptotic  theory  for  a  non-deformable  core 
approximation  o,  the  results  of  the  theory  for  a 
deformable  core  fo-  a  uniform  vortlclty 
distribution  o,  and  for  a  quadratic  distribution 
♦,  and  the  experimental  results  of  Widnall  et  al. 


In  Figure  10.  we 
number  n*  against  the 


plot  the  critical  wave 
self-induced  velocity  V, 


used  to  characterize  the  ring,  for  t*e  four  cases. 


On  the  same  pl.ne  we 
Widnall  et  al.  [15]  and 


reproduce  the  results  of 
Widnall  and  Tsai  [31]  for 


the  non-deformable  core  model,  the  deformable  core 
model  and  their  experimental  results.  The 
comparison  is  interesting  and  proves  our  early 
speculation  that  tne  numerical  thin  tube  model 
allows  small  core  deformation  since  the 
computational  results  are  closer  to  the 
experimental  data  than  those  of  the  analytical 
solution  of  the  non-deformable  core  model. 
However,  it  does  not  allow  enough  changes  within 
the  co-e  to  capture  higne-  order  radial  variations 
within  the  core  whicn  support  the  short  wave 
instability  that  is  observed  experimentally. 

So  far,  it  car  be  concluded  that  although  the 
results  of  the  tnin  tube  approximation  are  in 
agreement  with  the  analytical  theory  of  a  vortex 
ring  with  a  non-deformable  core,  the  model  is  not 
capable  of  describing  the  stability 
characteristics  of  a  vortex  ring  with  deformable 
finite  core,  as  observed  experimentally .  Using 
vortex  balls  allows,  however,  small  first  order 
deformations  in  the  vortlclty  core  of  the  ring,  as 
shown  in  Fig.  11,  which  move  the  predictions  of 
the  unstable  modes  closer  to  the  experimental 
values  than  the  analytical  theory  of  non- 
deformable  core  but  not  as  close  as  the  results  of 
the  more  elaborate  theory  of  deformable  cores. 
Thus  we  must  proceed  to  a  more  detailed 
description  of  the  vortlclty  core  of  the  ring 
using  several  vortex  balls  to  discretize  the 
vortlclty  within  the  core,  as  we  will  show  next. 
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Figure  11,  A  comparison  between  the  vortlclty 
distribution  within  the  vortex  ring  before 
deformation  and  after  deformation  of  amplitude  c/R 
•  0.1  and  ny  •  6  for  p/R  •  0.25,  computed  at  point 

of  zero  deformation. 


II. 5.  STABILITY  OF  A  VORTEX  TORUS 

To  make  the  distinction  between  the  two 
models  of  the  vortex  ring  clear,  this  ring  is 
called  a  vortex  torus.  In  this  representation, 
the  core  of  the  vortex  torus  is  discretized  into 
more  than  one  vortex  ball  so  that  6  <  0.  Thus, 
the  vortex  torus  is  formed  or  a  number  or  vortex 
rings  whose  cores  are  smaller  than  the  core  of  the 
torus.  The  initial  vortlclty  m,(0)  associated 

with  each  vortex  ball  is  computed  from  Eq.  (12)  by 
solving  the  corresponding  system  of  linear 
equations,  subject  to  the  condition  that  the  total 
circulation  is  the  same.  Since  the  torus  is 
uniform  in  the  azimuthal  direction,  it  suffices  to 
solve  a  number  of  equations  equal  to  the  number  of 
balls  used  across  one  cross-section  of  the  ring. 

In  the  results  presented  here,  nine  balls 
were  used  across  each  section  or  the  ring,  one  at 
the  center  and  eight  distributed  along  the 
circumference  of  a  circle  with  radius  p  •  O.T«o, 
arranged  at  *5*.  This  choice  for  the  initial 
location  of  the  centers  of  the  vortex  balls  is 
used  to  minimize  the  difference  between  the  total 
circulation  of  the  vortex  torus  and  the  sum  of  the 
circulation  of  the  vortex  balls.  The  core  radius 
of  each  ball  was  taken  as  «  -  0.8  0.  and  the 
distance  between  the  centers  or  neighboring 
elements  is  h  •  4/1.1.  Therefore,  the  number  of 
elements  used  along  the  circumference  of  the  torus 
depends  on  its  radius.  Figure  12  show  the  actual 
vortlclty  distribution  and  tne  numerical  vortlclty 
distribution  produced  by  the  vortex  balls.  The 
motion  of  these  balls  throughout  the  cross  section 
or  the  torus  allows  a  substantial  deformation  of 
its  core  at  different  sections.  Thus,  hlghe- 
order  radial  modes  associated  with  the  instability 
or  vortex  rings,  which  has  been  observed 
experimentally  and  analytically,  can  be  captured. 
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Figure  12.  The  actual  and  discretized  vorticity 
distribution  of  a  vortex  torus  using  nine  vortex 
balls  within  the  core. 


Computations  were  performed  for  four  tori 
with  o/R  -  0.15.  0.20,  0.25,  and  0.35.  In  all 

cases,  the  vorticity  distribution  Is  the  same 
third-order  Causslan  utilized  before.  The  initial 
amplitude  of  the  perturbation  c/R  -  0.02,  and  a 
number  n  of  slnewaves  were  fitted  along  tne  torus 
to  model  the  initial  perturbation.  The  time  step 
of  lnteg-ation  at  -  0.1,  and  the  computations  were 
performed  for  1000,  o-  2000  time  steps  to  observe 
the  growth  of  the  perturbation.  The  overall 
behavior  of  the  vortex  torus  was  the  same  in  all 
cases.  As  an  example,  results  for  o/R  -  0.2  are 
discussed  next  in  detail. 

Figure  >3  shows  two  views  of  the  torus  after 
1000  time  steps,  when  perturbed  with  n  •  8.  9.  10, 
and  11.  In  the  fi-st  three  cases,  the  core 
deforms,  as  seen  by  the  redlstrlbut ion  of  the 
individual  rings  wltnln  tne  torus,  and  the  waves 
spin  around  the  unpe-turbed  axis  of  the  ring. 
However,  the  perturbation  stays  bounded.  In  the 
last  case,  the  perturbation  grows  in  both  the 
radial  and  the  streamwlse  directions  causing 
substantial  non-uniform  deformation  around  the 
ring  axis,  and  the  generation  of  streamwlse 
vorticity.  The  amount  of  deformation  in  each  case 
is  seen  from  the  total  number  of  elements  used  at 
the  last  time  step.  In  the  first  three  cases,  the 
number  of  elements  remains  constant  at  N  •  1080 
for  1000  steps.  In  the  unstable  case,  the  number 
of  elements  grows  to  1262.  Since  from  Helmholtz 
theorem,  «(t)/«(0)  •  ai(t)/atf0),  where  atj  is  the 

summation  over  aij  for  all  the  vortex  elements, 

this  stretch  is  accompanied  by  intensification  of 
the  vorticity  within  the  ring  at  the  same  ratio  of 
stretch. 

Figure  1 A  shows  two  views  for  the  torus  in 
the  unstable  case  every  200  computational  time 
steps,  starting  at  t  •  0.  It  is  clear  that,  at 
the  unstable  mode,  waves  do  not  rotate  around  the 
a. is  of  the  ring  while  their  amplitudes  grow. 
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Figure  13.  The  form  of  the  vortex  torus,  e/R  • 
0.2,  after  1000  computational  time  steps  with  at  - 
0.1,  starting  with  a  perturbation  or  c/R  •  0.02. 
The  wavenumber  of  the  perturbation  is  n  •  8,  9, 
10,  and  n,  starting  from  the  top  plot. 
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Figure  1*.  The  growth  of  the  perturbation  on  * 
vortex  torus  with  e/R  •  0.2,  t/R  •  0.02  and  n  - 
n.  The  torus  is  discretized  into  9  rings  with 
120  vortex  balls  along  each  ring.  Results  are 
shown  every  200  time  steps  starting  with  t  •  0. 
The  plots  show  that  this  is  an  unstable  node  in 
which  the  amplitude  of  the  waves  grow  without 
rotation. 


similar  to  the  results  of  the  thin  tube 
app-oxlmatlon  and  to  the  results  of  Wldnall  and 
Tsai  [31].  Moreover,  the  core  vortlclty  Is 
redistributed  Into  a  number  of  sectors  equal  to 
the  number  of  waves.  The  outer  portion  of  each 
sector  stretches  forward  In  the  streaawlse 
direction  while  the  Inner  part  elongates 
backwards.  On  the  other  hand,  results  for  n  •  9, 
which  is  a  stable  node,  depicted  every  300  steps 
In  Fig.  15,  show  the  rotations  of  the  waves  as 
peaks  and  valleys  exchange  their  locations  while 
the  core  vortlclty  remains  uniform  In  the 
azimuthal  direction. 

The  average  anplltude  of  the  perturbation 
around  the  circumference  of  a  torus  with  o/R  - 
O.15  and  0.35  Is  shown  In  Figs.  16  and  17, 
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Figj-e  ’6.  The  growth  of  the  perturbation  with 
time  fo-  a  torus  with  o/R  •  0.2,  perturbed  with  n 
•  8 ,  9.  10,  and  11. 


respectively.  The  exponential  growth  of  snail 
perturbations  at  the  early  stages  is  seen  at  n  ■ 
11  and  6  for  the  two  eases,  respectively.  Figure 
18  summarizes  the  results  for  the  four  tori,  o/R  * 
0.15.  0.2.  0.25  and  0.35.  showing  a  very  good 

agreenent  with  the  experimental  results  of 
Widnall  and  Sullivan  [1*].  As  before,  the  value 

of  V  is  used  to  characterize  the  ring  In  order  to 
renove  any  confusion  about  the  definition  of  the 
core  and  the  vortlclty  distribution.  The 
analytical  results  for  a  vortex  torus  with  a  non- 
unlforn  vortlclty  distribution  within  the  core, 
the  nunerlcal  results,  and  the  experlnental  data 
are  In  close  agreement,  differences  can  be 
primarily  attributed  to  the  vortlclty  distribution 
within  the  core. 

The  form  of  the  unstable  torus  with  o/R  * 
0.35  at  n  •  6,  Is  shown  at  time  steps  1ROO-2OO0, 
every  200  steps  in  Fig.  19  (It  was  found  that  n  » 
7  Is  also  an  unstable  mode  for  this  torus).  It  Is 
Interesting  to  note  that  the  core  deformation  Is 
different  at  different  azimuthal  locations  and 
that  the  Inner  and  outer  radii  do  not  follow  the 
same  pattern  (Yule  [32]).  The  figures  Indicates 
that  the  Inner  and  outer  edges  of  the  vortlclty 
core  of  the  torus  may  move  in  anti -phase  and  that 
variations  at  frequencies  different  than  the 
perturbation  frequency  arise  at  late  times.  Thus, 
nigher  order  radial  modes  form  as  part  of  the 
Instability  of  vortex  rings,  In  accordance  with 
the  conclusion  or  the  analytical  theory  [15].  To 
quantify  these  frequencies,  we  study  the  energy 
spectra  of  two  tori.  In  Figs.  20  and  21,  the 
spectra  for  e/R  •  0.20  at  t  •  100  and  for 
o/R  •  0.35  at  t  •  200  are  shown.  In  the  stable 

modes,  only  the  perturbation  frequency  Is  present 
at  very  small  amplitude.  In  the  unstable  modes, 
nigher  harmonics  of  the  perturbation  frequency, 
which  had  zero  amplitudes  at  t  >  0  are  excited  at 
substantial  levels. 
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Figure  20.  The  frequency  spectra  for  the  energy  of 
a  torus  with  o/R  •  0.2  measured  at  t  ■  100,  for  n 
•  8**;n»9,»Oi  n  »  10  and  n  •  11  *0. 
In  the  unstable  node,  higher  harmonics  at  22  and 
33  a re  excited. 
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Figure  21.  The  frequency  spectra  for  the  energy  of 
a  torus  with  o/R  -  O.35,  measured  at  t  •  200,  for 
n  *  5  *  *  and  n  •  6  *  0. 
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Figure  19.  The  evolution  of  the  unstable  mode  of  a 
vortex  torus  with  e/R  •  0.35,  perturbed  at  t  •  0 
with  c/R  •  0.02  and  n  •  6.  The  torus  Is  shown  at 
time  steps  1*00-200,  every  200.  The  total  number 
of  vortex  elements  at  these  plots  Is  8 1 0  (as  at  t 
•  0),  822,  870  and  1032,  respectively. 


Results  of  this  section  Indicate  that  the 
numerical  simulation  of  vortex  rings  with  finite 
and  deformable  cores,  as  represented  by  a  number 
of  vortex  balls  distributed  within  the  ring  core 
and  along  its  axis,  can  accurately  be  used  to 
compute  the  wavenumber  of  the  unstable  modes  and 
their  growth.  The  deformation  of  the  ring  Into  a 
number  of  sectors  aligned  with  the  unstable 
standing  waves  resembles  strongly  the  experimental 
results.  The  generation  of  small  scales, 
accompanied  by  higher  frequencies  In  the  energy 
spectrum,  Is  due  to  vortex  stretch  and  leads  to 
the  well-known  turbulence  cascade.  Thus,  the 
azlsuthal  Instability  of  vortex  ring  may  be  one 
mechanism  of  transition  to  turbulence. 
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Ill  AXISYMMETRIC  SHEAR  LAYER 

The  ultimate  goal  of  this  work  is  to  study, 
using  three  dimensional  vortex  simulation,  the 
evolution  of  a  turbulent  axisymmetric  shear  layer. 
Analytical  and  experimental  studies  of  this 
configuration  Indicate  that  several  types  of 
instability  can  arise  and  influence  the 
development  of  the  flow  field  (e.g.  Yule  [32], 
Crow  and  Champagne  [33]  and  Mlchalke  and  Hermann 
[3*0.)  These  instabilities  Include  axisymmetric 
modes,  such  as  the  rollup  and  pairing  of  ring-like 
vortex  eddies,  and  the  Jet-preferred  mode,  as  well 
as  azimuthal  modes  such  as  the  type  which  was 
analyzed  in  the  previous  section.  The  interaction 
between  these  modes,  which  have  not  been  fully 
understood,  govern  the  flow,  and  in  particular, 
the  velocity  fluctuations,  entrainment,  and  mixing 
between  the  fluids  on  both  sides  of  the  layer. 
Chemical  reactions,  sound  generation,  combustion 
instability  can  be  strongly  affected  by  these 
interactions. 

In  this  section,  results  for  the  evolution  of 
an  axisymmetric  shear  layer,  subject  to  an 
axisymmetric  perturbation,  and  using  the  three- 
dimensional  vortex  scheme  developed  in  the 
previous  section,  are  presented.  The  computations 
are  restricted  to  a  periodically  excited  layer, 
thus  boundary  conditions  on  both  sides  of  the 
computational  domain,  l.e.  the  wavelength,  are 
satisfied.  This  is  accomplished  by  using  image 
vortices  on  both  sides  of  the  domain  and  computing 
their  field  by  summing  over  the  Induced  velocity 
of  these  Images.  The  first  image  of  each  vortex 
on  both  sides  must  be  considered  as  a  vortex  ball 
with  a  finite  core  radius.  Beyond  that,  images 
can  be  considered  as  vortex  points  with  zero 
cores.  The  effect  of  the  latter  can  be  arranged 
as  a  summation  ove'-  an  infinite  series  for  a  two- 
paramete-  function.  This  function  is  computed, 
using  a  large  number  of  terms  in  the  series,  and 
stored  as  a  two-dimensional  table  of  the  two 
parameters.  During  the  computations,  an 
interpolation  procedure  is  used  to  obtain  the 
value  of  the  function  at  Intermediate  points. 
Details  will  be  published  elsewhere. 

Results  are  obtained  for  an  axisymmetric 
shear  layer  with  the  following  parameters:  the 
thickness  of  the  vortlclty  layer  I/D  •  0.2, 
wavelength  of  the  perturbation  i/D  •  1.32,  and 
amplitude  of  perturbation  e/D  •  0.0k,  where  I  •  2 
o  and  a  is  the  standard  deviation  of  the  second- 
order  Gaussian  vorticity  distribution  within  the 
layer  and  D  is  the  mean  diameter  of  the  layer. 
The  layer  is  discretized  into  16  sections  in  the 
streamwlse  dl-ectlons  and  5  sections  in  the  cross¬ 
stream  direction,  resulting  in  80  vortex  rings. 
Each  ring  is  represented  by  50  vortex  balls  along 
its  axis.  The  vorticity  of  each  vortex  ball  was 
obtained  as  before  using  4/D  •  0.0825.  Figure  22 
shows  a  comparison  between  the  actual  and 
discretized  vorticity  distribution. 

Figure  23  shows  the  location  of  the  vortex 
elements  in  p-z  plane,  where  p  is  the  radial 
direction  and  z  is  the  streamwlse  direction,  and 
their  streamwlse  velocity  relative  to  the  mean 
flow.  In  Fig.  2k,  vortex  balls  which  were 


X 

Figure  22.  Actual  and  discretized  vorticity  within 
the  shear  layer.  Five  vortex  balls  are  used  at 
the  Indicated  location. 

neighbors  in  the  z-direction  at  time  t  •  0  are 
connected  using  cubic  spline  curves  to  show  the 
stretch  that  the  flow  experiences  while  vorticity 
is  during  the  rollup  phase.  The  vorticity  field 
and  the  material  lines  are  plotted  every  20 
computational  time  steps  starting  at  t  •  0. 
Although  five  layers  of  vortex  elements  were  used 
to  discretize  the  vorticity  in  the  radial 
direction,  the  plots  in  Fig.  2*  show  only  two 
layers,  the  central  layer  and  the  next  layer  to 
the  outside.  Plots  of  vortex  elements  locations 
in  the  radial  plane  and  the  p-z  plane  show  that 
the  elements  remain  on  perfect  circles  while  the 
radii  of  these  circles  Increase  or  decrease  as  the 
vorticity  layer  rolls  up. 

Results  in  Fig.  23  indicate  that  the  initial 
perturbation  causes  the  layer  to  rollup,  forming  a 
large  scale  ring-like  vortex  eddy.  As  time 
progresses,  more  of  the  vorticity  becomes 
concentrated  around  this  eddy,  and  more 
irrotatlonal  fluid  from  both  streams  is  entrained 
into  its  core.  Due  to  the  self-induced  velocity 
of  curved  vortex  lines,  the  eddy  moves  in  the 
streamwlse  direction.  However,  within  the 
duration  of  the  computation,  it  maintains  perfect 
azimuthal  symmetry.  Figure  2k  shows  that  the 
central  layer  experiences  the  strongest  stretch 
within  the  core  as  it  endures  several  turns  due  to 
secondary  instabilities,  while  the  "braids',  l.e. 
the  two  sleeves  connecting  neighboring  cores, 
become  thinner  due  to  the  strain  field  of  the 
cores. 

Since  the  layer  maintains  a  perfect 
axisymmetric  configuration  during  rollup,  one  can 
make  a  preliminary  conclusion  that  the  growth  of 
the  axisymmetric  modes  during  the  early  stages  of 
development  suppresses  the  azimuthal  instability 
modes  of  the  evolving  vortex-ring  eddy.  This  is 
widely  supported  by  the  analytical  linear  theory 
or  Mlchalke  and  Hermann  [35]  and  by  the 
experimental  results  of  Vandsburger  et  al.  [16], 
Roquemore  et  al.  [17]  and  Namazian  et  al.  [35]. 
The  analytical  study  shows  that  the  exponential 


Figure  23.  Rollup  of  an  axisyanetric  shear  layer 
due  to  an  axlsymmetrlc  perturbation.  Vortex  balls 
are  represented  by  circles  and  their  velocity 
relative  to  the  mean  velocity  Is  depleted  by  a 
line  vector  starting  at  the  center  of  the  circle. 
Results  are  plotted  every  20  time  steps  starting 
at  t  •  2.0.  The  solid  line  Indicates  the 
centerline  of  the  layer. 


Figure  2H.  Lines  connecting  neighboring  vortex 
balls  In  the  streanwise  direction.  The  figure 
shows  the  central  layer  and  the  next  layer  to  the 
outside.  The  solid  line  Indicates  the  centerline 
of  the  layer. 
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IV  CONCLUSIONS 


growth  rates  of  axisymmetric  modes  are  higher  than 
those  of  azimuthal  modes.  Moreover,  in  all  the 
reported  experimental  results,  Including  those  of 
Yule  [32]  and  Crow  and  Champagne  [33].  the  layer 
starts  6y  rolling  up  into  a  perfect  vortex  ring. 
At  later  stages,  the  rings  become  susceptible  to 
the  azimuthal  instability  discussed  in  the 
previous  section  and  the  flow  loses  the  azimuthal 
coherence. 

However,  this  Is  only  a  preliminary 
conclusion  since.  In  the  numerical  solution,  the 
flow  started  with  a  large  axlsymmetrlc 
perturbation  and  zero  azimuthal  perturbation.  The 
amplitudes  of  the  perturbations  were  selected  to 
model  most  experimental  situations,  where 
azimuthal  perturbations  are  Inhibited  at  the  onset 
of  the  layer  by  the  nozzle.  The  higher  growth 
rate  of  the  axisymmetric  mode  could  also  be  a 
property  of  the  linear  range,  as  shown  by  the 
analytical  results.  To  support  this  result  and  to 
study  the  effect  of  the  azimuthal  modes  on  the 
growth  of  the  axisymmetric  modes,  we  are  planning 
to  run  the  same  case  for  different  ratios  of  the 
axlsymmetrlc/azlmutnal  amplitudes.  Another  issue 
to  be  Investigated  is  at  what  stage  does  the 
azimuthal  modes  start  to  grow  and  what  effect  does 
the  strain  field  generated  during  pairing  have  on 
Its  development. 

Apart  from  the  displacement  of  the  ring  eddy 
in  the  streamwlse  direction  due  to  the  curvature 
of  tne  vortex  lines,  this  eddy  resembles  the  eddy 
that  forms  during  the  rollup  of  the  planar  shear 
layer  (Chomem  et  al.  [36]).  This  similarity  has 
bee-  observed  before  In  the  two-dimensional 
calculations  of  Davis  and  Moore  [37].  As  Fig.  2« 
indicates,  the  streamwlse  strain  field,  associated 
with  the  non-Iinea-  stages  of  rollup,  pulls  the 
vortex  elements  apart  so  that  the  distance  between 
the  neighbo-ing  centers  exceeds  h  by  a  large 
factor.  In  o-de-  to  maintain  the  -esolutlon  at 
late-  times,  the  vortieity  field  must  be 
redistributed  between  a  different  set  of  elements 
which  a-e  organized  so  that  they  can  accomodate 
this  strain  field,  as  suggested  by  Ghonlem  et  al. 
[36]  in  the  2D  vortex  element  method. 

The  stability  of  this  ring  eddy,  which  has  an 
elliptical  core,  to  azimuthal  perturbations  while 
it  is  moving  in  the  strain  field  of  Its 
neighboring  eddies  is  of  central  Importance  to  the 
later  stages  of  development  of  the  layer.  A 
numerical  study  of  this  problem  Is  unde-way. 


In  this  work,  we  have  embarked  on  the  task  of 
applying  three  dimensional  vortex  methods  for  the 
numerical  solution  of  the  unsteady  Navler-Stokes 
equations  at  high  Reynolds  number.  The  3D  vortex 
element  method  presented  here  combines  the 
advantages  of  the  vortex  filament  method  of 
Leonard  [8,n],  and  the  vortex  stick  method  of 
Chorln  [23,23]  in  maintaining  the  connectivity  of 
the  vortieity  field,  thus  satisfying  the 
solenoldallty  condition.  It  utilizes  the  results 
of  the  convergence  analysis  of  Beale  and  Majda 
[19,20]  in  selecting  the  core  of  the  vortex 
elements.  The  scheme  is  Lagranglan,  Is  capable  of 
capturing  the  effect  of  plain  strain  as  well  as 
the  vortex  stretching  along  vortex  lines  by 
changing  the  number  and  strength  of  the  vortex 
elements.  It  is  readily  extendable  to  flow  fields 
with  boundaries. 

Results  for  the  stability  of  a  vortex  ring 
with  a  finite  core,  which  forms  as  an  axisymmetric 
vortieity  layer  rolls  up,  show  very  good  agreement 
with  the  analytical  and  experimental  results.  The 
results  reveal:  (1)  the  breakdown  of  the  azimuthal 
coherence  of  the  ring  due  to  the  growth  of  radial 
perturbations  along  and  within  the  core;  (2)  the 
evolution  of  streamwlse  vortieity  In  the  non¬ 
linear  stages  of  Instability  In  the  form  of 
elongated  lobes  of  vortieity  along  wedges  within 
the  expanding  core;  and  (3)  the  development  of  an 
energy  cascade  to  small  scales  which  accompanies 
the  stretch  of  vortieity  during  the  non-llnear 
growth  of  Instability.  Similar  configurations 
were  captured  In  experimental  studies  on  vortex 
rings  and  later  stages  of  turbulent  Jets. 

The  scheme  was  used  to  Investigate  the 
Initial  stages  of  transition  to  turbulence  in  an 
excited  axlsynmetrlc  mixing  layer,  and  the  results 
showed  good  agreement  with  the  recent  qualitative 
results  of  Vandsburger  et  al.  [16]  and  Roquemore 
et  al.  [IT].  Quantitative  study  will  be  performed 
to  Investigate  the  Interaction  between  the 
axisymmetric  and  the  azimuthal  Instability  modes 
and  their  effect  on  the  development  of  the  flow. 


ACKNOWLEDGEMENT 

This  work  Is  supported  by  the  Air  Force 
Office  of  Scientific  Research  Grant  AFOSR  Sk-0356, 
The  National  Science  Foundation  Grant  CPE-8A0A81 
and  the  Department  of  Energy  Contract  DE-ACOk- 
66AL16310.  Computer  support  Is  provided  by  grants 
from  NSF  and  the  John  von  Neumann  Computer  Center, 


16 


REFERENCES 


Van  Dyke,  M.,  An  Album  of  Flute  Motion, 
Parabolic,  176  p.,  1982. 

Lugt,  H.J.,  Vortex  Flow  in  Nature  and 
Technology.  Wiley,  297  p. ,  '983. 

Ho,  C.-M.,  and  Huerre,  P.,  "Perturbed  free 
•hear  layers",  Ann.  Rev.  Fluid  Mech..  i_6,  pp. 
365-*2*  1 198V) . 

Robert,  R.A.  and  Roshko,  A.,  "Effects  of 
periodic  forcing  on  mixing  in  turbulent  shear 
layers  and  wakes,  AIAA-85-0570. 

Chonlem,  A.F.  and  Ng,  K.K.,  "Numerical  study 
of  a  forced  shear  layer",  Phys.  Fluids.  1987, 
in  press,  AIAA-86-0370. 

Chorln,  A.J.,  "Numerical  study  of  slightly 
viscous  flow,"  J.  Fluid  Mech,.  57.  pp.  *23“ 
**2  (1973). 

Saffman,  P.G.  and  Merlon,  D.I.,  "Difficulties 
with  three-dimensional  weak  solutions  for 
lnviscld  incompressible  flow",  Phys.  Fluids, 
29.  pp.  2373*2375  (1986). 

Leonard,  A.,  "Computing  three  dimensional 
incompressible  flows  with  vortex  elements", 
Ann.  Rev.  Fluid  Hecn..  17.  pp.  525*559 

TWrsT. 

Ashurst,  W.T.,  "Numerical  simulation  of 
turbulent  mixing  layers  via  vortex  dynamics, 
Proc.  the  2ist  Sym.  Turbulent  Shear  flow,  p. 
*02.  Sprlnger-Verlag,  1979. 

Ghonlem,  A.F.  and  Sethlan,  J.A.  "Effect  of 
Reynolds  number  on  the  structure  Of 
recirculating  flow",  AIAA  J.  1987,  in  pres*. 
AIAA-85-01V6. 

Leonard,  A.,  "Vortex  methods  for  flow 
simulation”,  J.  Comput.  Phys..  37,  pp.  289* 
335  (i960). 

Saffman.  P.G.,  and  Baker,  G.R.,  "Vortex 
interactions",  Ann.  Rev.  Fluid  Mech. .  11,  pp. 
95-122  (1979). 

Saffman,  P.G.,  "The  velocity  of  viscous 
vortex  rings".  Stud.  Appl.  Math..  *9.  oo. 
371-380  (1970).  *“ 

Wldnall,  S.E.  and  Sullivan,  "On  the  stability 
of  vortex  rings",  Proc.  R.  Soc.  Lond..  A  332, 
PP.  335-353  0973). 


16.  Vandsburger,  H.  Lewis,  G.,  Seltxman,  J.M., 

Allen,  M.G.,  Bowman,  C.T.  and  Hanson,  R.K., 
"Flame-flow  structure  In  an  acoustically 

driven  Jet  flame".  Western  States  Section/The 
Combustion  Institute.  19B6  Fall  meeting? 

Paper  14-1 4. 

17.  Roquemore,  W.M.,  Bradley,  R.P.,  Jackson, 
T.A.,  Klztrnls,  S.W.,  Gross,  L.P.,  Switzer, 
G.L.,  Trump,  D.D.,  Sarka,  B.,  Ballal,  D.R., 
Llghtman,  A.J.,  Taney,  P.P.  and  Chen,  T.H., 
"Development  of  Laser  diagnostics  for 

combustion  research",  Central _  States 

Section/The  Combustion  Institute.  19^  Spring 
Meeting. 

18.  Batchelor,  G.K.,  An  Introduction  to  Fluid 
Dynamics.  Cambridge  University  Press,  6'5  P? 

Tirr 

19.  Beale,  J.T.  and  Majda,  A.,  "Vortex  methods  I: 
Convergence  In  three  dimensions",  Hath. 
Comput. .  £9,  pp.  1-27  (1982). 

20.  Beale,  J.T.  and  Majda,  "Vortex  methods  II: 

Higher  order  accuracy  In  two  and  three 
dimensions,  Math.  Comput.,  pp.  29-52 

(1982). 

21.  Truasdall,  C.,  The  Kinematics  of  Vortlclty, 
Indiana  Univarsity  Prass,  232  p.,  195V. 

22.  Anderaon,  C.  and  Greengard,  C.,  "On  vortex 
methods”,  SIAM  J.  Numer.  Anal.,  22,  pp.  *»1 3— 
MO  (1985). 

23.  Chorln,  A.J.,  "Vortex  models  and  boundary 
layer  instability  SIAM  J.  Scl.  State. 
Comput..  1_,  pp.  1-21  (1980). 

2V.  Chorln,  A.J.,  "Estimates  of  tntermlttency , 
spectra  and  blow-up  In  developed  turbulence, 
Commun.  Pure  Appl.  Math..  jV ,  pp.  853*866 
(i$8iT 

25.  Chorln,  A.J.,  "The  evolution  of  a  turbulent 
vortex",  Commun.  Math  Phys.,  83.  PP.  517-535 
(1982). 

26.  Leonard,  A.,  "Numerical  simulation  of 

interacting,  three  dimensional  vortex 
rilaments",  Proc.  of  the  *th  Int.  Conf. 
Numer.  Meth.  Fluid  Dyn..  pp.  2*5-50. 

Sprlnger-Verlag  (1975). 

27.  Ashurst,  W.T .  and  Melburg,  E.,  "Three 

dimensional  shear  layers  via  vortex 
dynamics",  Sandta  National  Laboratory  report, 
SAND85-8777.  1985. 


Wldnall,  S.E. ,  Bliss,  D.B.  and  Tsai,  C.-Y., 
"The  Instability  of  short  waves  on  a  vortex 
ring",  J.  Fluid  Mech..  60,  pp.  35-*7  ('97*). 


28.  Greengard,  C.,  Three  Dimensional  Vortex 
Methods.  Ph.D.  thesis,  University  57 
California,  Berkaley,  198*. 


29.  Beale,  J.T.  and  Majda,  "Vortex  methods  II: 
Higher  order  accurate  vortex  methods  with 
explicit  velocity  kernels”.  J.  Comput.  Phys.. 
58,  pp.  >83-208  t 1985) . 

30.  Hidnall,  S.E.,  "The  structure  and  dynamics  of 
vortex  filaments,"  Ann.  Rev.  Fluid  Mech. .  pp. 
141-165  ;'976). 

31.  Wldnall,  S.E.  and  Tsai,  C.-Y.,  "The 

Instability  of  the  thin  vortex  ring  of 
constant  vorticity",  Philos.  Trans.  R.  Soc. 

'  Lond.  Ser.  A.  .,287.  pp.  273*305  0977)  . 

32.  Yule,  A .J..  "Large-scale  structures  in  the 
mixing  layer  of  a  round  Jet".  J .  Fluid  Mech. . 
89,  pp.  413-432  (1978). 

33.  Crow,  S.C.  and  Cnampagne,  F.H.,  "Orderly 
structure  in  Jet  turbulence",  J.  Fluid  Mech. . 
48,  pp.  547-591  (1971). 

34.  Michalxe,  A.  and  Hermann,  0.,  "On  the 

lnviscid  instability  of  a  circular  Jet  with 
external  flow",  J.  Fluid  Mech. ,  1 14.  pp.  343- 
359  (1982). 

35.  Namazian,  M.,  Kelly,  J .  Schaefer,  R . , 

Johnston,  S.  and  Long,  M.,  "Flow  and  mixing 
structures  in  bluff-body  stabilized  gas 

flames",  1986  International  Gas  Research 

Conference,  Toronto,  Canada. 


36.  Gnoniem,  A.F.,  HeldarineJad,  C.  and  Krlsnnan, 
A.,  "Vortex  element  simulation  of  the  rollup 
and  mixing  in  a  thermally-stratl f ied  shear 
layer",  submitted  for  publication,  J.  Comput. 
Phys. .  1987. 

3".  Davis.  R.W.  and  Moore,  E.F.,  "A  numerical 
study  of  vortex  merging  in  mixing  layers", 
3nys.  Fluids.  28,  pp.  1626-1635  ('985). 


t’j  7  . 

i <•  - :  ’  "  ' 


18 


,  ?t  °6T  y 


